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Abstract 

We address multiscale elliptic problems with random coefficients that are a perturbation of multiscale determin- 
istic problems. Our approach consists in taking benefit of the perturbative context to suitably modify the classical 
Finite Element basis into a deterministic multiscale Finite Element basis. The latter essentially shares the same ap- 
proximation properties as a multiscale Finite Element basis directly generated on the random problem. The specific 
reference method that we use is the Multiscale Finite Element Method. Using numerical experiments, we demon- 
strate the efficiency of our approach and the computational speed-up with respect to a more standard approach. We 
provide a complete analysis of the approach, extending that available for the deterministic setting. 

C<| 1 Overview of our approach and results 
in 

^-H The Multiscale Finite Element Method (henceforth abbreviated as MsFEM) is a popular numerical approach for 
.multiscale problems (see [371 IMl ISOl [321 1311 131 13S1 123 UHl HI])- It consists in a Galerkin approximation of the original 

T-H 'problem over a finite dimensional space generated by basis functions that are specifically adapted to the problem under 
I consideration. 

I ' This approach is popular for a tvi^ofold reason. First, its use is not restricted to multiscale problems that converge 
^ to a homogenized problem in the limit of vanishing ratio between the small scale and the macroscopic scale. It may be 
applied to much more general situations. Second, when the problem does converge to a homogenization problem, the 
MsFEM approach is meant to approximate the solution of the problem with the small scale e at its actual small value 
and not "only" in the asymptotic regime £ — >■ 0, which is the regime addressed by homogenization theory. 
To fix the ideas, consider the problem of finding u'^ solving 

- div [A^Vu^] = / in I?, = on dV, (1) 

on a bounded domain D C M'', with / G L^{'D), and where is a uniformly bounded, coercive matrix that varies at 
scale e. A standard Finite Element Method (FEM) would require a space discretization of the domain at the scale e 
in order to capture the oscillations of at scale e. This is prohibitively expensive. The MsFEM aims at accurately 
approximating using a limited number of degrees of freedom. It does not require the matrix A"^ to be periodic 
(namely A'^{x) — Aper{x/e) for a fixed periodic matrix Ap^r) or stationary. 

We now briefly describe the approach and present the aim of this article. Starting from a coarse mesh Th with a 
standard (say Pi) Finite Element basis set of functions {4'^}._^, generating the associated space 

Vh :== span(0°,i = 1, • • • ,L), 

we first numerically build the MsFEM basis functions 4>f . Several definitions of these basis functions have been proposed 
in the literature (yielding different numerical methods), and we detail this in the sequel (see e.g. pII)) - (ITT]) - (IT^ '). For 
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the moment, it is sufficient to know that, to each which varies at the macroscopic scale, is associated a function 
with variations at the scale e. In practice, 0? is numerically computed (in fact, pre- computed), using the specificities of 
the problem addressed. These highly oscillatory functions generate the finite dimensional space 

Wh :=span(0^,i = I,-- - ,L). 

Note that Wh and Vh share the same dimension. 

We next define the MsFEM solution um using a Galerkin approximation of ([T]) on Wh, instead oiVh- Again, details 
will be given below. The MsFEM solution um provided by the approach reads 

L 

1=1 

for some coefficients {{UM)i}f^i- Of course, these coefficients depend on e, but this dependency is kept implicit in the 
sequel. 

We now turn our attention to the stochastic problem 

- dW [A%-, uj)Vu%-, u)] ^ f inV, u^-,uj) ^ on dV, (2) 

and a typical quantity of interest E [u'^(a;, •)], which is traditionally approximated using a Monte Carlo method. In- 
troducing a set of A4 realizations of the stochastic matrix {^^'™}i<m<Ai' ^ direct, naive application of the MsFEM 
paradigm would consist in first computing for each realization m the stochastic MsFEM basis functions 0^'™(x,a;), 
next performing a Galerkin approximation of ([2]) using this MsFEM basis set to compute {w^(a;, and 
eventually approximating E [u^{x, •)] by 

M 

m=l 

Such an approach is unpractical because of the prohibitively expensive computational load. 

To reduce the computational cost and make the MsFEM approach practical in such a stochastic context, a natural 
idea we investigate in this article is to consider a less generic setting, for which a dedicated, more computationally 
affordable approach, can be designed. One possibility is to consider matrices {x, w) = A'^ (x) +B{x, uj) in ^ that are 
not highly oscillatory in their stochastic part. In such cases, dedicated approaches have been proposed, we refer to [3S] 
for more details. Another approach is to reduce the number of Monte-Carlo simulations used for the computation of 
the multiscale basis functions. In [l][26], the authors assume that their coefficient can be written as a Karhunen-Loeve 
type expansion, and apply a collocation method to a priori choose some sparse realizations for which they compute 
the multiscale basis functions. 

In this article, we consider one of the many alternate variants of problem We suppose that A^{x,uj) is highly 
oscillatory in both its deterministic and stochastic components, but that it is a perturbation of a deterministic matrix. 
More precisely, we assume that 

A'{x,Lu) = A'^{x,uj) = Al{x) +r]Al{x,uj), (3) 

where Aq is a deterministic matrix and ?7 is a small deterministic parameter. This model may be well suited for 
heterogeneous materials (or, more generally, media) that, although not periodic, are not fully stochastic, in the sense 
that they may be considered as a perturbation of a deterministic material. We call this setting the weakly stochastic 
setting. Note that many practical situations, involving actual materials or media, can be considered, at a good level of 
approximation, as perturbations of a deterministic (often periodic) setting (see e.g. [41]). 

In a series of recent works (see [Tll[Tni[2S] and [SUZIIH]; see also [S] for a unified presentation), we have considered 
such a setting, in the context of homogenization theory (the matrix A^{x,uj) in Q-© reads A'^{x,uj) = A^(a:/e,w) for 
a stationary matrix A^(a;,w), which is, in a sense to be made precise, a perturbation of a periodic matrix). We have 
shown there that, in such a case, the workload for computing the homogenized solution is significantly lighter than 
for generic stochastic homogenization, and actually comparable to the workload for periodic homogenization. We will 
show in the sequel that the MsFEM can be adapted to this weakly stochastic setting, providing an approximation of 
the solution to (H))-®, for fixed e, at a much smaller computational cost than the direct approach. 
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The main idea of our proposed approach is to compute a set of deterministic MsFEM basis functions using Aq , the 
deterministic part of Af^ in the expansion ([3]), and then to perform Monte Carlo realizations at the macroscale level 
using a set of Ai realizations of the random matrix w)}^^^^^^^ (see Section [2] for a detailed presentation). 

Note that, for each of these realizations, we solve the original problem, with the complete matrix A^^, and not only its 
deterministic part. Only the basis set is taken deterministic. By construction, the approach provides an approximation 

L 

i=l 

of u^^{x,uj), where the basis functions are deterministic. These basis functions are computed only once, hence the 
cost to compute {^"(a;, w)}^^^^^^ is much smaller than the cost to compute {u']^{x,uj)}^^^^j^. This is especially 
true if ^ has to be solved for many right-hand sides /. We expect that this approximation us is as accurate as um 
for small rj. We show below that this is indeed the case, when A^^ is a perturbation of Ag (see Section [3] for numerical 
tests). 

We would like to note that the MsFEM is not the only multiscale technique based on finite elements. The bottom 
line of our approach, consisting of generating suitable multiscale functions for the discretization of a weakly stochastic 
problem, using for this purpose the deterministic reference problem, can in principle be applied to other multiscale 
techniques. Another popular technique is the HMM method [27[ [25] . for which our approach could in principle be 
easily adapted. 

In the numerical tests reported on in Section [3j we compare, in the norm, (the exact solution to ([2]) with 
the matrix A^ = given by (|3])) with us (the solution provided by our approach) and um (the solution provided by 
the ideal, expensive approach). The quantity — um\\h^{v) somewhat represents the best possible accuracy that we 
can achieve, in the sense that our approach inherits the limitations of the MsFEM approach. We thus cannot expect 
our approximation us to be more accurate than um- We can only hope to compute an approximation of comparable 
quality with a much reduced workload. The numerical results we obtain confirm that, for small ry in ([3]), the quantity 
\\us — uli\\m{v) is of the same order of magnitude as \\um — ul^\\m['D)i although, we repeat it, the computational cost 
to compute us is much smaller than that to compute um- 

We next derive error bounds for our approach in Sectional We recall that, in the deterministic setting, a classical 
context for proving convergence of the MsFEM approach is the case when, in the reference problem (IT|), the matrix 

reads A^{x) — Aper ^— ^ for a fixed periodic matrix Ape^- Likewise, to be able to perform our theoretical analysis in 

the stochastic setting, we assume in Section 2] that A'^^ix^ui) — Ay^ ^— ,a;^ for a fixed stationary random matrix A^. 

The problem ©-([SI) then admits a homogenized limit when e vanishes. 

Our proof follows the same lines as that in the deterministic setting, which we now briefly review (see the introduction 
of Section 2] for more details on the structure of the proof) . The MsFEM is a Galerkin approximation, that we assume 
momentarily, for the sake of clarity, to be a conforming approximation (this is indeed the case when, for defining the 
highly oscillatory basis functions , oversampling is not used) . The error is then estimated using the Cea lemma: 

ll^t" - '"m||_h-i(X)) < C* inf Wu" - VhWmm), 
vueWh 

where is the solution to the reference deterministic highly oscillatory problem ([T]), um is the MsFEM solution and 
C is a constant independent of e and h. Taking advantage of the homogenization setting, we introduce the two-scale 
expansion 



i=l 



of zi^, where u* is the homogenized solution, w^. is the periodic corrector associated to G M'', and diU* denotes the 
partial derivative — — . We next write 

OXi 



um\\hi(v) <C [Wu" -v''\\hi{v)+ inf Wv" - Vh\\m{v}] 
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The first term in tlie above riglit-liand side is estimated using standard homogenization results on the rate of convergence 
of — u^. To estimate the second term, one considers a suitably chosen element Vh S W/i, for which — Vh\\H'^{'D) 
can be directly bounded. 

Following the same strategy in our stochastic setting, we estimate the distance between the solution uf^ to the 
reference stochastic problem (H))-® and the weakly stochastic MsFEM solution Ms as 



,uj) - us{-,uj)\\h^v) < C - v^^{-,uj)\\Hi{v) + ^i^i^^ - Vh\\H^v)^ 



We observe that a key ingredient for the proof is the rate of convergence of the difference between the reference 
solution ufj and its two-scale expansion w^. Such a result is classical in periodic homogenization, but, to the best of 
our knowledge, open in the general stationary case (in dimensions higher than one). One only knows that wf^ — 
vanishes (in some appropriate norm) when e — > 0. However, in the particular case when is only weakly stochastic, it 
is possible to obtain such a result, as we have shown in [?^. Hence, exploiting the specificity of our weakly stochastic 
setting, we are able to obtain (see our main result, Theorem [TU] and estimate (p5|) ): 
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< 



C (^V^ + h+^+?] (^^y^^ \n{N{h)) + 77 + r]^C{?])^ 



where || • ||^i is a broken norm, C is a constant independent of £, h and 77, C is a bounded function as 77 goes to 
0, and N{h) is the number of elements in the mesh (roughly of order /i^** in dimension d). As is often the case in 
the deterministic setting, we use here (both for our numerical tests and in the analysis) the oversampling technique. 
Consequently, the basis functions do not belong to Hq{T>), hence the use of a broken norm in the above 
estimate. As we point out below, when 7; = in ([S]) , our approach reduces to the standard deterministic MsFEM (with 
oversampling), and the above estimates then agree with those proved in |32) . 

This article is organized as follows. First, in Section[2l we describe the MsFEM approach. For consistency, we begin 
by the deterministic setting in Section 12. 1[ and point out there that the direct adaptation to the general stochastic 
setting yields a prohibitively expensive approach. The adaptation of the approach to the weakly stochastic setting 
is described in Section 12.21 We next turn to numerical simulations, in Section [31 Some procedures to efficiently 
implement the approach are first described in Section [XTl We next consider a one-dimensional test (see Section 15^ . 
which is useful for several reasons. First, it allows to calibrate some numerical parameters, such as the number M of 
independent realizations when estimating the exact expectation by an empirical mean. Second, we assess the accuracy 
of our approach with respect to the magnitude of 77. We demonstrate there that 77 does not have to be extremely small 
for our method to be very efficient. For instance, on the test case considered in Section [X^ we show that our approach 
is as accurate as the expensive, direct approach as soon as 77 is such that 

ml 



is equal to or smaller than 0.1, 

ixn) 

where Oq is the deterministic component of the diffusion coefficient and 770^ is the stochastic component (see 
expansion ([3])). Lastly, we also assess the accuracy of our approach with respect to the presence of frequencies in the 
random coefficient that are not taken into account in the MsFEM basis set. We next turn to two test cases in 



dimension two, where we observe that our approach behaves as well as in the one-dimensional case (see Section [ 
In particular, in Section r3. 3. 21 we successfully address a classical test-case of the literature. 

Section |4] is devoted to the analysis of the approach, in the homogenization setting. Our main result. Theorem IIOI 
is presented in Section HTTl and proved in Section The proofs of some technical results are collected in Appendix YK\ 
In addition, we specifically consider the one dimensional case in Section [ 



2 MsFEM-type approaches 

For consistency and the convenience of the reader, we present in this section the MsFEM approach to solve the original 
elliptic problem ([1]). For clarity, we begin by presenting the approach in a deterministic setting. The reader familiar 
with the MsFEM may easily skip this section and directly proceed to Section 12.21 where we present our approach in a 
weakly stochastic setting. 
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2.1 Description in a classical deterministic setting 

Let u"^ e H^{T>) be the solution to ([1]), where the matrix A"^ £ {L°° {'D)Y^'^ satisfies the standard coercivity condition: 
there exists two constants a_|_ > a_ > such that 

Ve, VxeP, VCeM'', a^\£_\^ <£_'^A'{x)^ and \\A'\\l^(^v) < a+. 

Note that the MsFEM approach is not restricted to the periodic setting. We therefore do not assume that A'^{x) = 
Aper{x/e) for a fixed periodic matrix ^per- 

The MsFEM approach consists in performing a variational approximation of ([1]) where the basis functions are 
precomputed and encode the fast oscillations present in ([l}. In the sequel we argue on the following formulation, 
equivalent to ([T]): 

Find e H^{V) such that, for any v G H^iV), Aeiu", v) = b{v), (4) 

where 

x)Vu(x) dx and b{v) — / f{x) v{x) dx. 
Jv Jv 

The MsFEM is a three-step approach: 

1. introduce a standard discretization of the domain V using a coarse mesh as compared to the small scale oscillations 
of A^. 

2. for each element K of the coarse mesh, compute the basis function (j)^' as the solution of an elliptic equation 
posed in K (see e.g. (Il0l)-(in])-(IT2l) below). 

3. solve the Galerkin approximation of (|H), for the set of basis functions defined at Step 2. 

The advantage of the approach is that, for the same accuracy of the approximation as that provided by a standard 
FEM, the macroscale mesh can be chosen sufficiently coarse so that the resulting discretized problem has a limited 
number of degrees of freedom, and may thus be computationally solved inexpensively. This is observed in practice |37| , 
and proven by a theoretical analysis (see [3H1I22]) when the problem (U) admits a homogenized limit. See also 50! and 
references therein. 

To further illustrate this fact, we reproduce here a simple one-dimensional analysis we borrow from A. Lozinsky 
(see [m Chap. 6] and [IZ]). This analysis explains remarkably well the interest of the approach, and, in contrast 
to [55J 121], is not restricted to a homogenization setting. Consider the one-dimensional domain V = (0, 1) and the 
reference problem 

Cu = f, it(0) = 'u(l) = 0, 

for the operator Cu — (i/u')', where / e L^(0,1) and ly € L°°(0, 1) with i'{x) > Vmin > almost everywhere on 
(0, 1). The function may have oscillations at a small scale. The associated weak formulation reads 

Find u e Hq{0, 1) such that, for any v e Hq{0, 1), a{u, v) = b{v), (5) 

with 



a{u,v) ^ / h'{x)u' {x)v' (x) dx and b{v) = / f{x)v{x)dx. 
Jo Jo 

We now introduce the nodes = a;o < xi < ■■■ < xl — 1 that define the elements Ki — [xi--i,Xi]. Let h = 
max \xi — Xi-i \ be the mesh size. The multiscale finite element space 

Wh = {vh e C°(0, 1) such that Cv,, = on each K^} , (6) 

defined using the operator C, is adapted to the problem under study. We next proceed with a Galerkin approximation 
of (0 using the space Wh- 

Find Uh e Wh such that, for any Vh G Wh, a{uh,Vh) = b{vh). 
The solution Uh then satisfies 

\\u- Uh\\E < —J==\\f\\L^OA) (7) 
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where || • Hb = \J a(-, •) is the energy norm. The proof of this estimate goes as fohows. By definition of u and Uh^ we 
have a(u — it/i,w/i) = for any Vh G W?i. Hence, Uh is the orthogonal projection of u on W/i according to the scalar 
product a(-, •). Since || • \e is the norm associated to that scalar product, we have 

\u-Uh\E^ inf u/illfi. (8) 

Choose Vh to be the finite element interpolant of u, which is defined by Vh(xi) — u{xi) for any i — 0, 1, . . . , L, and 
consider the interpolation error e = u ~ Vh- On each element Ki^ we have, precisely because the space Wh is defined 
as (ISl), 

Ce — — (i/e')' — f with e{xi^i) — e{xi) — 0. 
We multiply by e, integrate by part and obtain 

i^{x)\e'{x)\'^ dx ^ f{x)e{x)dx <\\f\\L2(^K,)\\e\\L^Ki)- (9) 

1 Jxi-i 

Since e vanishes on the boundary of Ki, the Poincare inequality with the best constant {xi — Xi-i)/TT yields 
Ml^ko < ^ iWWmK,} < ^ — / i'ix)\e'ix)f dx 



By substitution in ([5]), we obtain 

|2 ^ "• II -f II 2 



Kx)|e'(a:)|2dx<-5^||/||i.(^^). 

Summing over the elements and using ([5]) yields ([7]). Using again that v is bounded from below, we deduce from ([T]) 
that ^ 

\W-Uh\\H^0.1) < 7; II/IIl2(0,1)' 

where C-p is the Poincare constant of the domain V = (0, 1). As pointed out in [33J Chap. 6], the interest of the above 
estimate (or of estimate ([7])) lies in the fact that the constant in the right-hand side only depends on v through Vmim 
and remains the same even if v oscillates at a small scale. In contrast, for a standard finite element method, the error 
is also proportional to /i, but with a constant that depends on the norm of the exact solution u. With a standard 
finite element space Wh, we indeed classically deduce by Cea's lemma that 

II II , II'^IIl~(0,1) ■ r II II lkllL=°(0,l) II n II 

where C-p is the Poincare constant of the domain T) = (0, 1), and RhU is the projection of u on Wh according to the 
scalar product. We thus obtain that 

II II , ^, lkl|L~(0,l) II „2 II 

where C is independent from the functions v and u. If v oscillates at a small scale (e.g. v{x) = vixje) for a fixed 
function i/), the iP' norm of u may be large (of the order of e""'^). A FEM approach then requires h to be smaller than 
e to reach a good accuracy. 

We conclude this illustration by noting that such a general analysis of the MsFEM approach is not available in 
dimension d > 2. The analysis presented in |38l I32j . which is performed without any restriction on the dimension, 
additionally assumes that the matrix in Q reads A^(x) = Aper{x/e) for a fixed periodic matrix Aper- 

We now describe the MsFEM in a multidimensional setting. 
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Definition of the coarse mesh For simplicity (see Remark [T] below) , we consider a classical Pi discretization of 
the domain T). We denote by Th the corresponding mesh, with L nodes. Let (/)°, z = 1, • • • , L, be the basis functions. 
We introduce the finite element space 

Vh ■■= span(0°,i = 1, • • • ,L), 

and define the restriction 

^0,K ,01 

of these functions in each element K. 

Remark 1. We refer to for a presentation of a MsFEM method that uses P2 macroscale basis functions. 



Definition of the MsFEM basis Several definitions of the MsFEM basis functions have been proposed in the 
literature (see e.g. [37l |38l [32l [3l [30l [S^). They all follow the same pattern but they give rise to various methods. We 
present in the following the particular method that we have implemented. It makes use of the oversampling technique 
introduced in [37] and developed in 

For any element K, we consider a domain S D K (see Figure [T]), obtained from K by an honiothetic transformation 
of center the centroid of K, and of ratio larger than 1. 




xf 



Figure 1: Definition of S (in 2D for clarity) 

Let denote the coordinate of the vertex j of the domain S. For any vertex i of S, we introduce the affine function 

Xi'^ (defined on S) that satisfies the condition Xi'^i^f) — ^ij for J- Let xt^ G H^{S) be the unique solution to the 
problem 

(10) 



A^x)Vxt^ix) -0 inS, x-'^=X^^ on 38, 



which, in practice, is numerically solved e.g. using a finite element method with a mesh size adapted to the small 
scale £. We then define the local basis functions 



d+1 

,e,K \ ^ e,S 



K 



(11) 



as linear combinations of the restrictions of xl'^ on K, with aij chosen such that 



d+1 



VI < i,j < d + 1, 



,0,K/ K 



(12) 



where xf denotes the coordinate of the jth vertex of the element K. Note that the condition (|T2|) is enforced on the 



function (j)^'^ , and not on 4>1''^. The coefiicients aij are consequently independent from e. As (ji^''^ and x 
both affine on K, condition (IT^ implies that 



0,K 



0,S 



K 



d+1 



Vl<i<d+1, VxeK, .^-^"^(x) -^ay■x°'^(^)■ 
J=l 



(13) 
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We next introduce the functions defined on 2? by '/>f Ik = for elements K. 

Note that the problems (jlOp . indexed by S, are all independent from one another. They may be solved in parallel. 

Macroscale problem Wc now introduce the finite dimensional space 

Wh :=span(</)^, i = l,--- ,L), 

and proceed with the approximation 

Find um G VV/i such that, for any v S Wh, A^{um,v) ~ b(v), (14) 

of dU, where 

A^,{u,v)^y] [ {\7v{x)fA''{x)Vu{x)dx and b{v) = [ f{x)v{x)dx. 

Observe that i/)^ has jumps across the edges of the triangulation (due to the use of the oversampling technique), hence 
Wh <t H^{'D), thus the broken integral used to define A^{u,v). On the other hand, since Wh C L^{'D), the linear 
form b is well defined for v € Wh- The formulation (jl4p is a non- conforming Galerkin approximation of This 
brings additional error terms in the error estimation (see Lemma |12l in Section U]). On another note, remark that the 
dimension of Wh is equal to L. The formulation (jl4p hence requires solving a linear system with only a limited number 
of degrees of freedom. 

We are now in position to substantiate our claim in the introduction, where we briefly mentioned that, in the 
stochastic setting, a direct application of the MsFEM to approximate the solution to ^ is unpractical. It would indeed 
lead to compute, for each realization of A'^(x,ll)), first a basis set and second a macroscale solution. This approach has 
been briefly examined theoretically in [21. It is prohibitively expensive. We therefore turn to an alternate approach. 

2.2 A weakly stochastic setting 

We now restrict the general setting and propose a dedicated, practical MsFEM type approach. Following up on 
previous works (see [HI [SI [Ml EI] ) and as announced in (jH]), we assume here that the random matrix A^{x,u}) in ^ 
is a perturbation of a deterministic matrix, in the sense that 

A%x,u;) = A^ix,u;) = A'^ix) + A^x , uj) , (15) 

where 77 g R is a small deterministic parameter, Aq and are bounded matrices, and is coercive, uniformly in e. 
We also assume that the matrix A'^ itself satisfies the coercivity and boundedness assumptions, uniformly in rj and e 
(we refer to jH [71 [H] and |15[ for other perturbative settings) . 

The principle of the proposed approach is to compute the MsFEM basis set of functions with the deterministic part 
Aq of the matrix Af^, and next to perform Monte-Carlo realizations for the macroscale problem (|^- (fT5l) . where we keep 
the exact matrix A'^ (and not only its deterministic part). Following the approach sketched in Section l2.ll we first 
solve l|10p . with A'^(x) = Aq(x), and build the deterministic finite dimensional space 

Wh ■■= span(0f, I = 1, • • • ,L) 

following (fTT|) - ([T^ . We next proceed with a standard Galerkin approximation of ([^- ([T51) using Wh- For each m €E 
{1, • • • ,M}, we consider a realization A'^'"^{-,uj) and compute u™(-,a;) G Wh such that 

VweW/,, V/ {\/v{x)fAf^'"'{x,uj)\7u'^{x,uj)dx^ [ f{x)v{x)dx. (16) 
Jk Jv 

Since the MsFEM basis functions are only computed once (rather than for each realization of A^^{x,lu)), a large 
computational gain is expected, and obtained, in comparison to the direct approach described above. 
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3 Numerical simulations 



This section is devoted to the many numerical simulations we have performed. We first discuss some implementation 
details. Next, we numerically estimate the performance of our approach on various test cases, and assess its sensitivity 
with respect to the magnitude of r/. We consider in Section 13.21 a test case in dimension one. In Section 13.31 we next 
study two test cases in dimension two. We also study how the presence in Af (the random component of the matrix 
A"^) of high frequencies that are not present in the deterministic component Aq, and that are thus not encoded in the 
highly oscillatory basis functions, affects the accuracy of our approach. 

Let ufj be the reference solution to Q-® obtained using a finite element method with a mesh size adapted to the 
small scale e, us be the approximation given by our approach (described in Section [2.2p and um be the approximation 
given by the direct approach (in which the MsFEM basis set is recomputed for each realization A^'™'{x, w), as explained 
at the end of Section l2.ll) . Our goal is to compare the error us — of our numerical approximation with the error 
uai — W^ri of direct and expensive approach. When 77 is small, we expect the approximation us to be essentially as 
accurate as the approximation um, and we show below that this is indeed the case. 

In the sequel, we assess the accuracy using the estimators 

. X ^ f \\ui-U2\\LHv) \ , , . ^ [ \\ui~U2\\hi \ 
eL2{ui,U2) =E[ — - — ^ and ej^i (ui, U2) = E — - — , (17) 

V \\u2\\l^v) ) \ \\u2\\hI ) 

where u\ and U2 are the solutions obtained with any two different methods, and 

1/2 

:= I V Mmm 1 (18) 



II«IIhi(k) 



is the broken li^ norm. The expectation is in turn computed using a Monte-Carlo method. Considering M realizations 

||mi(-, w) — U2(-, w)||jyl 

{_Xraiy>y\ \<^.m<M '^'^ ^ random variable, e.g. X(bj) = — -j- ^, we compute the empirical mean [Im and 

- - ll"2(-,t^)||//i 



the empirical standard deviation om as 

^ A/ , M 



A/ A/ 



m— 1 m—l 

As a classical consequence of the Central Limit Theorem, the following estimate is commonly employed: 

\E{X)-^,MiX)\<lM 

It provides a practical evaluation of E(X) from the knowledge of ^m{X) and auiX). The numerical parameters 
have been determined by an empirical study of convergence. For instance, for the reference solution, we choose the 

mesh size h such that the quantity — is smaller than 0.03 %, thereby formally admitting that the 

11"^' \\m{-D) 

approximation has converged in h. The MsFEM parameters are determined likewise. 

All the computations have been performed using FreeFem++ 33 , with the MPI tools. 



3.1 Implementation details 

In the deterministic version of the MsFEM, the same matrix A'^ appears in the definition (jlOp of the basis functions 
and in the macroscale variational formulation (jl4p . This can be used to expedite the computation of the stiffness 
matrix associated with p^ . In our approach, described in Section the matrix that appears in the definition of the 
basis functions is Ag, whereas the macroscale variational problem involves = Aq + rjA\. An additional numerical 
computation is thus needed. 

To solve (ITSl) . we need to compute, for each element K and each realization A^'™(a;,aj), the integrals 

^ir^= I (v#^(a;))^A^^^(a;,c.)V(/.5'^(a:)dx, (20) 
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where cj)^'^ are deterministic functions. We recall that Af^{x, lu) — AQ{x) + riAl {x, w) (see (llSp ). To allow for an efficient 
evaluation of (|20l) . we assume henceforth that is of the form 



AUx,u:) = J2 iQ+fe (7) Xkiu;) B^ix), 



(21) 



where Q = (0,1)'*, where {Xk)f.^id are scalar random variables, and for any k <E Z"*, x 1— >■ B^{x) G W^'^'^ are some 
deterministic functions. We comment on this assumption in Remark [2] below. The important consequence of (j2ip is 
that we can write the integral (1201) as a linear combination of deterministic integrals over cells of size e, with random 
coefficients. To simplify the notation, we assume that the spatial dimension is d = 2. We define 



P 



q 



(m Vj_ Vk_ 

\ e e e 



1, 



• (Vi V3 Vk^ 

— ,— ,— 

V £ e e / 

and likewise, we define the integers / and m (see Fig. [5]). We can then write (QUI) as 

q—l m—1 



K 



(V0^'^(x)) A^^"\x,u:)Vc^)^{x) dx = K^lf" + 



a=p j3=l 



where 



0,K 



1,K 



(v</)r^(a:))^A^(a:)V0^^'^(a:)dx, 

k(x) (v#^(x))'^i?,^2;)V0^^'^(x) 



K 



dx. 



(22) 

(23) 
(24) 




Figure 2: To practically compute the integral (I^Ul) . we write that each element K (here in dimension ci = 2) is a subset 
of a quadrangle (here [le, me] x [pe, qe]) composed of cells of size e'^. 

We thus compute once the deterministic integrals and (IMl) . Next, for each realization of Af^, we evaluate 

the stiffness matrix elements /C^^-'"(a;) using the right hand side of ([2^ . No numerical quadrature is needed. As a 
consequence of (|2ip . most of the work for assembling the stiffness matrix is only performed once, independently of the 
number of Monte Carlo realizations. This significantly contributes to the gain in term of computational cost. 
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Remark 2. Assumption (j2ip is quite general, and already covers many interesting cases in practice. As explained 
above, the point in (j2ip is that A\ is a direct product ( or here, a sum of direct products ) of a function depending on x 
with a random variable that only depends on uj. Otherwise stated, A\{x,uj) depends linearly, in an explicit way, of lo. 
A similar assumption is made when applying reduced basis methods JJ^ to a problem of the type 

Find u\ such that, for any v, a{u\,v;X) = b{v), (25) 

where a(-,-; A) is a bilinear form parameterized by A. Assume this problem has been solved for some values {Ai}[^-^ 
of the parameter, yielding the functions {u\-Y-^^. Under the assumption that a(-,-;A) — ao(-, •) + Aai(-,-) (namely, 
a(-, •; A) depends linearly on X), one can precompute the stiffness matrix elements ao{u\^ , u\- ) and ai(u\. , u\.) for any 
1 < < ^- This allows to next perform a very efficient Galerkin approximation of the problem (PS)) (for any X) on 
the space Span{ux. ,i — 1, - ■ ■ , I). 



3.2 One-dimensional test-case 



The purpose of this section is threefold. We first calibrate the number M of realizations considered for the Monte-Carlo 
method for the two-dimensional numerical experiments that we consider in the sequel. We next investigate how the 
accuracy of our approach depends on rj and on the presence of frequencies in the random coefficient that are not 
taken into account in the MsFEM basis set functions. The low computational costs that we face in this one-dimensional 
situation allow us to test our approach more comprehensively than in the two-dimensional test-cases described below. 

Let (^fe)/,g2 denote a sequence of independent, identically distributed scalar random variables uniformly distributed 
in [0, 1]. We consider the random coefficient 



a^(a;, ^ik,k+i] (-) ( 5 + 50 sin^ ) -|- T]Xk{uj) k sin^ 

which is a particular example of the expansion (jl5p with 

al{x,uj) 



(^nx 

e 



al{x) = 5 + 50sin^ ^'^"^ 
and that satisfies the structural assumption (PT|) 

K.{k,0 -- 



Y'^{kM+i] (-) Xk{uj) K sin^ f^^j , 
fcez ^ \ ^ / 

We set e — 0.025 and choose k such that the quantity 

af(-,tj) 



— SupEss, 



(26) 



L~(X)) 



has the same value JC — I ioi the three different values oi ( = {1, 3, 7} we consider below. We analytically compute the 



reference function u^, solution to 



d 

dx 



du'i 



1 in (0,1), ii^(0,c^) = 1*^(1,0.) =0, 



as well as the MsFEM basis functions for both approaches. Let um and us be the approximation of by the two 
MsFEM approaches described above, where the coarse mesh size is h = 1/30. 

We first calibrate the number of independent realizations to accurately approximate the exact expectation in (|17p 
by the empirical mean (jl9p . To this aim, we present on Fig. [3] the mean and the confidence interval computed using 
(|19l) for an increasing number M of realizations (we compute up to 1000 independent realizations). We check that this 
indicator reaches a plateau for M > 30, and thus converges fast. On this example, considering 30 realizations is hence 
sufficient to accurately compute the error ([T7|. Based on this observation, we will only consider AI — 30 realizations 
in the two dimensional examples of Section 13.31 



Remark 3. There is no reason to think that the calibration of our parameters that we perform in the one- dimensional 
situation provides an adequate adaptation of these parameters for the higher dimensional setting. We however see no 
other manner to proceed and the approach has indeed provided us with good results. 

Note also that the MsFEM approach is much more accurate in the one- dimensional setting than in the two- 
dimensional setting (compare Tables{^ and\^ with Tables\^ and \10\ below). This is due to the specificity of the 
one dimensional setting. However, one- dimensional examples remain relevant for e.g. assessing how the MsFEM 
accuracy depends on rj. 
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I 

6- I 



I 

4 - I 
I 

2 - I 
I 

0- ' 

-2 1 ' ' ' ' ' ' ' ' ' ' 

5 10 15 20 25 30 35 40 45 50 

M 

Figure 3: Convergence of the indicator eHi{uM,u'^) (see PT|) ). for 77 = 1, C = 1 and k = 55. For each value of M, we 
plot the empirical mean along with its confidence interval, computed from the first M realizations. We only plot the 
results for the first 50 realizations. 

We now check how the accuracy of our approach depends on rj. In Tables [1] [31 [31 01 [S] and [SI we report the estimators 
(|17p. along with their confidence intervals, for various choices of (k, () that all correspond to IC{k, Q = 1. For rj < 0.1, 
we observe that \\us — u^\\ and \\um — are of the same order of magnitude, and are both larger than \\um ~ us\\ 
(both in and broken norms). We thus obtain the same accuracy with the direct and the weak stochastic MsFEM 
approaches, whereas the weak stochastic MsFEM is computationally (much) less expensive. For 77 = 1, as expected, 
the accuracy of the approximation us deteriorates. The accuracy of um is independent of 77. 

Remark 4. In Section^ we estimate in the (broken) norm the error between the reference solution u,^ and the 
weak MsFEM solution us- For information, we also include in Tables[Ii\S[ the numerical comparison in the norm. 

We now turn to a different question. In the example considered here, some frequencies present in af do not appear 
in Oq, and are thus not captured in the highly oscillatory basis functions 4>f. We now show that our approach can still 
handle this case, provided the amplitude of these modes remains small. 

We first consider the case when the amplitude k associated to the frequency ( is kept constant, and compare the 
performance of our approach in the case C — 1 and = 3. In the latter case, a relevant high frequency is not taken 
into account in the basis set functions. Comparing Tables [1] and [H (corresponding to C = 1) with Tables [3 and [51 
(corresponding to = 3) for a given value of 77, we see that the accuracy of our approach deteriorates. This is not 
unexpected, of course. Otherwise stated, to achieve a given accuracy (say an error of 15 % in the broken norm), 
we need to take smaller values of 77 (namely 77 < 0.01) when C = 3 than when C = 1 (in which case 77 — 0.1 is already a 
sufficiently small value). 

We now run the comparison differently. As we increase the gap between the frequency present in af and that 
present in ag (i.e., as we increase C), we simultaneously decrease the amplitude k of that mode. In practice, we do this 
by keeping constant the parameter /C(k, C) defined by (|26[). Then the accuracy of our approach remains constant, and 
is independent of C. See indeed the numerical results of Tables [T][6l that all correspond to the choice IC{k,() = 1, for 
three different values of We observe that, at fixed 77, errors are comparable, and independent of the value of (k, C). 

In conclusion, the accuracy of our approach depends both on the amplitude k and the value ^ of the high frequency 
not taken into account in the MsFEM basis set functions. If ( and k are scaled so that /C(k, () remains constant (which 
implies that k decreases if C increases), then the accuracy of our approach remains constant. 
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Table 1: g^(0, 1) error (fTT]) (in %) for k = 55 and C ^ 1 



V 




ei/i(M5,u^) 




1 

0.1 
0.01 


0.14644 ±0.00036 
0.16001 ±0.00006 
0.16258 ±0.00000 


2.62550 ± 0.02696 
0.15021 ± 0.00051 
0.10837 ± 0.00002 


2.44359 ± 0.02696 
0.07036 ± 0.00044 
0.04825 ± 0.00025 


Table 2: H^O, 1) error ^ (in %) for k = 14.38 and C = 3 


V 


e/fi(wAf,w^) 


e^i(w5,uQ 


em{us,UM) 


1 

0.1 
0.01 


0.18269 ±0.00030 
0.16529 ±0.00003 
0.16314 ±0.00000 


2.38950 ± 0.02277 
0.14959 ± 0.00055 
0.10840 ± 0.00000 


2.23869 ± 0.02230 
0.08082 ± 0.00041 
0.04954 ±0.00001 


Table 3: H^O, 1) error ^ (in %) for n = 8.39 and C = 7 






em ("s,<,) 


em{us,UM) 


1 

0.1 
0.01 


0.17436 ±0.00026 
0.16465 ±0.00004 
0.16308 ±0.00000 


2, 34495 ±0.02105 
0.15748 ±0.00067 
0.10846 ±0.00000 


2, 27358 ±0.02089 
0.09803 ± 0.00053 
0.05054 ± 0.00001 



Table 4: L^(0, 1) error p7)) (in %) for k = 55 and C = 1 









6^2(775, um) 


1 

0.1 
0.01 


0.00018 ±0.00000 
0.00018 ±0.00000 
0.00018 ±0.00000 


0.07286 ± 0.00317 
0.00045 ± 0.00002 
0.00015 ± 0.00000 


0.06861 ± 0.00306 
0.00024 ±0.00001 
0.00002 ± 0.00000 


Table 5: L^{Q, 1) error ^ (in %) for k = 14.38 and C = 3 






eL2(ws,u^) 


6^2(7*5, TiM) 


1 

0.1 
0.01 


0.00019 ±0.00000 
0.00018 ±0.00000 
0.00018 ±0.00000 


0.06658 ± 0.00270 
0.00036 ± 0.00001 
0.00015 ± 0.00000 


0.06238 ± 0.00261 
0.00019 ±0.00001 
0.00002 ± 0.00000 


Table 6: L^{0, 1) error ^ (in %) for n = 8.39 and C = 7 


'/ 




6^2(7/5, u^) 


6^2(775, TiAf) 


1 

0.1 
0.01 


0.00018 ±0.00000 
0.00018 ±0.00000 
0.00018 ±0.00000 


0.08903 ± 0.00310 
0.00037 ± 0.00002 
0.00015 ± 0.00000 


0.08410 ± 0.00261 
0.00016 ±0.00000 
0.00003 ± 0.00000 



Table 7: g^(0, 1) error ^ (in %) for k ^ 55 and ( = 3 



V 






eHi{us,UM) 


1 

0.1 
0.01 


0.21826 ±0.00073 
0.17142 ±0.00013 
0.16383 ±0.00001 


12.30047 ±0.10647 
0.59293 ± 0.00519 
0.11448 ±0.00014 


12.01694 ±0.10617 
0.49523 ±0.00489 
0.05247 ± 0.00007 


Table 8: L^{0, 1) error 1^ (in %) for k = 55 and C = 3 


V 


e£2(77Af,W^) 


eL2(ws,w^) 


6^2 {us,um) 


1 

0.1 
0.01 


0.00022 ±0.00000 
0.00019 ±0.00000 
0.00018 ±0.00000 


1.53780 ± 0.03878 
0.00503 ± 0.00027 
0.00018 ± 0.00000 


1.51837 ± 0.00385 
0.00406 ± 0.00024 
0.00005 ± 0.00000 
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3.3 Two-dimensional test-cases 

We now test our approach on two-dimensional test cases. Using the first test case, we show, similarly to the one- 
dimensional situation, that the weak stochastic MsFEM yields accurate results, provided the parameter rj is sufficiently 
small, and provided that the amplitude associated to frequencies present in Af^ but not encoded in the deterministic 
basis functions is small (see Section [3.3.ip . Next, in Section [3.3.21 we consider a test case similar to a classical benchmark 
test case of the literature. We again observe that our approach is efficient. For both cases, we show that the parameter 
r] does not need to be extremely small for our approach to be highly competitive. 

3.3.1 A multi- frequency case 

In line with what we observed in the one-dimensional case, we show here that the weak stochastic MsFEM provides 
interesting results even in the case when not all the frequencies present in are captured in the deterministic basis 
functions, provided their amplitude is not too large. To this aim, we consider the following numerical example. 

Let (-^fe,0(fe denote a sequence of independent, identically distributed scalar random variables uniformly dis- 
tributed in the interval [0, 1]. We consider the random matrix 

A^,j{x,y,uj) = al{x,y) Ids + Tyaf (x, y, cj) Idz, 

with 



ao{x,y) = 
a\{x,y,uj) 




Again, this choice is a particular example of the expansion (fT5|) satisfying the structural assumption (j2ip . We consider 
two different values of C, namely C = 1 ^-nd C = 3. As in the previous test case, the frequency ^ is not present in the 
deterministic part of Af^, and thus not encoded in the basis functions. In line with what we observed in Section 13.21 
we choose the amplitude k associated to that frequency such that the quantity (I26p has the same value K, = 1 for both 
values of C,. We compute solution to 

-div Vu^(-,cj)] = 1 in P, M^(-,w)==0 ondV, 

on the domain V = (0, 1)^ with e = 0.025. Let um and us be its approximation by the two MsFEM approaches 
described above. The numerical parameters for the computation are again determined using an empirical study of 
convergence. We use for the reference solution a fine mesh of size hf = e/40. The MsFEM basis functions are 
computed in each element K using a mesh of size /im — e/80. The oversampling parameter (i.e. the scale ratio of the 
homothetic transformation between K and S, see Fig. [Ij is equal to 3. The coarse mesh size is ft. = 1/30. In view 
of the results of Section 13.21 we consider M = 30 independent realizations, which will prove to again be sufficient to 
obtain accurate results. 

In Tables [9] and [To] (Tables [11] and [12] respectively), we report the estimator (fT7|) . along with its confidence interval, 
for the broken H^{'D) norm and for the i^(2?) norm, respectively. The results obtained here confirm our observations 
in the one-dimensional setting (Section 13. 2[) : 

• for given ^ and k, we observe that, when 77 is sufficiently small (here, rj < 0.1), the alternative approach provides 
a solution us that is an approximation of as accurate as um, for a much smaller computational cost (as the 
MsFEM basis set has only been computed once rather than for each independent realization of A^^). 

• our approach yields accurate results even if the frequency C is not encoded in the basis functions (/>f , provided the 
associated amplitude k is scaled accordingly. Figures in Table [H] (respectively Table [TT|) are very close to those of 
Table [To] (respectively Table [T2|) . This confirms that the error made by the weak stochastic MsFEM seems to be 
independent of k and ^, provided these two parameters are scaled so that JC{k, (^) remains constant. If ^ becomes 
different than 1, the frequency present in Oq, then the amplitude k associated to the frequency ^ has to decrease 
to keep /C(k, Q (and thus the accuracy of us) constant. 

These observations again demonstrate the efficiency of the approach. 
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Table 9: H\V) error 1^ (in %) for k ^ 73.61 and C = 1 



V 


e^i(MAf,<,) 


e^i(Ms,"^) 


e^i {us,um) 


1 

0.1 
0.01 


7.8437 ±0.1350 
6.8053 ±0.0165 
6.7338 ±0.0017 


19.8818 ±0.4123 
7.3868 ± 0.0276 
6.9795 ±0.0016 


18.8662 ±0.4216 
3.1528 ±0.0517 
1.8763 ±0.0013 


Table 10: H^{V) error ^ (in %) for k = 10 and C = 3 


V 






em {us,um) 


1 

0.1 
0.01 


6.7224 ±0.0368 
6.7154 ±0.0044 
6.1725 ±0.0004 


12.7292 ±0.2172 
7.1069 ±0.0128 
6.9770 ±0.0010 


10.8128 ± 0.2442 
2.2925 ± 0.0206 
1.8504 ±0.0003 


Table 11: L^{V) error ^ (in %) for n = 73.61 and C = 1 


V 






(us, Mm) 


1 

0.1 
0.01 


1.4355 ± 0.0795 
1.0630 ±0.0108 
1.0211 ±0.0011 


4.1649 ±0.1652 
1.1369 ±0.0075 
1.1512 ±0.0007 


2.8468 ±0.1694 
0.1441 ±0.0354 
0.1351 ±0.0014 


Table 12: L^{V) error (H?]) (in %) for k = 10 and C = 3 


V 


ei2(uAf,uQ 




6^2 (us, Mm) 


1 

0.1 
0.01 


1.0744 ±0.0127 
1.0226 ±0.0015 
1.0170 ±0.0001 


1.8433 ±0.0582 
1.1249 ±0.0038 
1.1551 ±0.0004 


0.8426 ±0.0832 
0.1147 ±0.0073 
0.1427 ±0.0003 



Remark 5. In TablesW ilSl we observe that the size of the confidence interval is much smaller than the distance between 
two different errors. This a posteriori validates the choice of the number M of Monte Carlo realizations according to the 
calibration we performed in the one- dimensional setting. In the two-dimensional setting studied here, we observe that 
considering M = 30 realizations is again sufficient. The same conclusion holds for results presented in Tables rOHTSI 
below. 

3.3.2 A classical test case 

We consider in this section a test case similar to a classical test case of the literature (see e.g. [371 [Ml UHl 132] )■ Let 
(^fe.j)(fc i)ez2 denote a sequence of independent, identically distributed scalar random variables uniformly distributed 
in the interval [0, 1]. We consider the random matrix 

AF/ N , f^\-, fy\ /2 + Psin(27ra:/£) 2 + sin(27rj//e) \ , ^^ ^ , 

^ni-^y^-^= E /(M+U y (f ) ( 2 + Psin(2.,/e) + 2 + Psin(2../e) j + ^^^-'^'^ 

with P ~ 1.8 and e — 0.025. We compute the reference solution u^ and its two approximations um and us with the 
same numerical parameters as in Section [3.3. II 

In Tables [T51 and [Til we report the estimator (|17p . along with its confidence interval, for the broken H^{T>) norm 
and for the L'^i'D) norm, respectively. We again see that, when rj is sufficiently small, us is an approximation of the 
reference solution m^ as accurate as mm- In Tables [T5l and [T6l we report on the accuracy of us, for more values of ry. 
Assuming that the accuracy of mm does not depend on 77 (which is consistent with the results reported in Tables [T31 
and [T^ . we see that our approach is as accurate as the direct, expensive MsFEM approach, as soon as 77 < 0.1 (if we 
use the broken norm to assess accuracy) and 77 < 0.25 (if we rather use the norm). The parameter 77 hence does 
not need to be extremely small for our approach to be highly competitive. 
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Table 13: H^(V) error ([TT]) (in %) 



V 


eHi{uM,ul) 


eHi(ws,u^) 


e^i {us,um) 


1 

0.1 
0.01 


8.1154±0.1913 
7.1664 ±0.0199 
7.1453 ±0.0020 


17.3678 ±0.7784 
7.0524 ±0.0705 
7.2837 ±0.0067 


15.5113 ± 0.8689 
2.5638 ±0.1006 
1.3882 ± 0.0020 


Table 14: L'^(V) error ^ (in %) 


V 




6^2(^5, U^) 


eL^{us,UM) 


1 

0.1 

0.01 


0.5620 ±0.0803 
0.5354 ±0.0160 
0.5347 ±0.0012 


1.6855 ±0.4860 
0.5688 ±0.0630 
0.6192 ±0.0054 


1.4739 ±0.5048 
0.1984 ±0.0712 
0.1072 ±0.0032 



Table 15: H^{V) error ^ (in %) 



Table 16: L^{V) error (d?]) (in %) 



'? 




V 


6^2(^5, U^) 


1 


17.3678 ±0.7784 


1 


1.6855 ±0.4860 


0.5 


15.9578 ±0.3461 


0.5 


1.0246 ±0.4414 


0.25 


10.6130 ±0.1591 


0.25 


0.5291 ±0.2285 


0.1 


7.0524 ±0.0705 


0.1 


0.5688 ±0.0630 


0.01 


7.2837 ± 0.0067 


0.01 


0.6192 ±0.0054 



4 Analysis 

This section is devoted to the analysis of the approach introduced in Section [2?2l and to the derivation of error bounds. 
As is often the case for the MsFEM (see e.g. [32]), we perform the analysis in a setting where the problem ©-([HI) that 
we consider admits a homogenized limit as e vanishes (although, we repeat it, the approach is used in practice for more 
general cases). The structure of our proof is similar to that for the deterministic setting, which we now overview (we 
refer to [32. for all the details). 

In the case when the oversampling technique is not used, the MsFEM is a conforming Galerkin approximation, and 
the error is estimated using the Cea lemma: 

IW" - umWh^ < C inf \\u^~Vh\\H^, 
vheWh 

where is the solution to the reference deterministic highly oscillatory problem ([1]), um is the MsFEM solution, and 
the constant C is independent from e and h. In the case when the oversampling technique is used, the MsFEM is a 
non-conforming Galerkin method. The error is then bounded from above by the sum of the best approximation error 
(the right-hand side of the above estimate) and the non-conforming error (that we do not detail here): 



umWh 



1 < C 



inf \\u'^ — t'/i||/fi + non-conforming error 



Note that, in the non-conforming case, the MsFEM solution um does not belong to H^, and one should write the above 
estimate with a broken norm rather than the norm. For the sake of clarity, we ignore this distinction in this 
preliminary discussion. 

Taking advantage of the homogenization setting, we introduce the two-scale expansion 



d 

=1 
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of u'^, where u* is the homogenized solution, ui^. is the periodic corrector associated to Ci £ W^, and diU* denotes the 
partial derivative — — . We next write 



dxi 



inf 

vheWh 



Vhlln^ + non-conforming error 



The first term in the right-hand side is estimated using standard homogenization results. To estimate the second term, 
one considers a suitably chosen element Vh G W/i, for which — w/iH^/i can be estimated directly. The main idea is 
that the highly oscillating part of can be well approached by an element in Wh, since, by construction, the highly 
oscillatory basis functions are defined by a problem similar to the corrector problem, and thus encode the same highly 
oscillatory behavior as that present in the correctors Wg. . We are thus left with approximating the slowly varying 
components of , for which standard FEM estimates are used. Lastly, we again use the fact that our problem admits 
a homogenized limit to estimate the third term, i.e. the non-conforming error. 



In the sequel, we follow the same strategy in our stochastic setting. We hence first write (see ([M]) below) that 



\\u^{;Lo)-usi;uj)\\m <C 



inf 



Vh{-,tL>)\\H^ + non-conforming error 



(27) 



where wf^ is the solution to the reference stochastic problem ©-(O and C is a deterministic constant independent from 
£, h and ij (note that, in (|64p . we use a broken norm rather than the norm; as pointed out above, this is due 
to the fact that our approach is a non- conforming Galerkin approximation; we ignore this distinction in the current 
discussion). To estimate the best approximation error (the first term in the right-hand side of ((27)) above), we use the 
triangle inequality, and write (see (l84)) below) that 



inf ||u^(-,w) 

Vh&Wh 



Vh{-,i^)\\H^ < IK,(-,'^) 



W)|ki 



(28) 



where is the two-scale expansion of the solution uf^ truncated at order e^. A first difficulty owes to the fact that, in 
the general stochastic setting, no estimate is known on ||M^(-,a;) — ti^(-, w)||jLfi. One only knows that its expectation 
vanishes when e — > 0. However, in the present article, we consider a weakly stochastic case. In that setting, we have 
derived such a convergence rate type result in 05] , and we can thus bound the first term of pSjl (see Section 14.1.21 
below for more details). The second term, inf ||u^(-,a;) — w/JIhI; of (128L is estimated using an explicit construction 

of a suitable Vh (see similarly to the deterministic setting. We again use there our specific weakly stochastic 

setting. Lastly, the non-conforming error (the second term in the right-hand side of (j27p above) is estimated following 
arguments similar to those of the deterministic case, using that our problem admits a homogenized limit and is weakly 
stochastic. 



This section is organized as follows. The error estimation is presented in Section 14.11 We first recall in Sec- 
tion mTTT] the formulation of the homogenized problem, and some results specific to the weakly stochastic case. Next, 
in Section I4.1.2[ we establish an error bound between the reference solution wf^ and its two-scale expansion (see 
Theorem [7]), which allows to bound the first term in the right-hand side of ([55)1 . Our main result. Theorem [TUl is given 
in Section [4.1.3l and proved in Section [4.21 The proof essentially consists in explicitly building a function Vh G such 
that the second term of ([28]) can be directly estimated. It also makes use of several technical results (Lemmas [TS] [Ml 
and 1161 below) to bound the non-conforming error, i.e. the second term in the right hand side of (j27p . The proof of 
these technical results is postponed until Appendix[Xl Last, in Section [4.31 we specifically consider the one dimensional 
case. 



Before proceeding further, we recall the setting of stochastic homogenization we work with. The reader familiar 
with this theory may directly proceed to Section [4.11 Let (r2,J^, P) be a probability space. For a random variable 
X G L^{U,,(]F), we denote by E(X) = j^-^X{ijj)dF{Lxj) its expectation value. We assume that the group (Z**, +) acts on 
ri. We denote by (Tk)kez,'^ this action, and assume that it preserves the measure P, i.e. 

Vfc e Z'^, VA e J", P{TkA) = P{A). 

We assume that r is ergodic, that is, 

e J", (Vfc e Z'^, TkA ^ A) ^ {¥{A) or 1). 
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We define the following notion of stationarity: any F S Ll^^ (]R'^,i^(r2)) is said to be stationary if 

yk e Z'^, F{x + k,uj) = F[x, TkUj) almost everywhere, almost surely. (29) 

Note that we have chosen to present the theory in a discrete stationary setting, which is more appropriate for our 
specific purpose, which is to consider a setting close to periodic homogenization. Random homogenization is more 
often presented in the continuous stationary setting. This is only a matter of small modifications. We refer to the 
bibliography for the latter. 

For the sake of analysis, we assume in this section that the matrix in ([l])-© reads A'^^{x,uj) ~ Ajj I—, 

where the random matrix Ajj is stationary in the sense of (1^^ . The problem ^ now reads 



e ' 



Ar. 



^(-,cj)v<,(-,w)] =/ inP, u^(-,cj) = ondV, (30) 



where A^{-,uj) G {L°°{W^)Y'"^ satisfies the standard coercivity and boundedness conditions: there exists two constants 
a+ > a_ > such that 

Vry, V^eM'', a_ < ^^(x, w)^ • ^ a.e. on M'^, a.s. and P^(-, w)||icc(Rd) < a+ a.s. (31) 

Due to the stationarity assumption on A^f, the problem (j30p admits a homogenized limit when e — 0. Note that, to 
the best of our knowledge, all analyses of the MsFEM approach in the deterministic setting that have been proposed 

in the literature are performed under a similar assumption (the matrix A^ in ([1} is assumed to read A^{x) — Ap^r ^— ^ 

for a fixed periodic matrix Ap^r, see e.g. [551 15^ ). 

In addition, in line with ([3]) and (ITS]) , we assume that A,, is of the form 

An{x,uj) ^ Aperix) +r] Ai{x,uj), (32) 

where S M is small parameter (we henceforth assume that \ri\ < 1), Ap^r is a symmetric bounded Q-periodic matrix 
{Q = [0, l]'') satisfying the ellipticity condition almost everywhere on R'', and Ai is a symmetric bounded stationary 
matrix: |Ai(x,a;)| < C almost everywhere in R"*, almost surely. Since r/ is small, our problem is weakly stochastic. 
In line with (I21L we furthermore assume that Ai is of the form 

Ai{x,Uj)^ ^ lQ+k{x)Xk{uj) Bper{x), (33) 

where (^^(tj))^.^^^ is a sequence of i.i.d. scalar random variables such that 

3C, Vfc e Z'^, \Xk{uj)\ < C almost surely, 

and Bper E (i°°(R''))'^^'' is a Q-periodic matrix. Besides being used in Theorem [7] below, this assumption is also 
used in the proof of Lemma I16[ to recognize that some quantity (namely, (I126P below) is a normalized sum of i.i.d. 
variables, on which we can use Central Limit Theorem arguments. As mentioned in Section 13.11 above, the form p3p 
is not essential. The point in (I33p is that Ai is a sum of direct products of a function depending on x with a random 
variable only depending on lo. Assumptions alternative to (j33l) could be made, that still satisfy this framework. 
Finally, we assume that 

Aper is Holder continuous, (34) 
Bper is Holder continuous. (35) 

We use these assumptions to obtain a rate of convergence of the two-scale expansion of (see [12 ^-nd Theorem [7] 
below), and hence control the first term in the right-hand side of (|28p . Such assumptions are standard when proving 
convergence rates of two-scale expansions (see e.g. [40l p. 28]). In turn, to control the second term in (j28l) and the 
non-conforming error (the second term in (j27p ). we do not need Bp^r to be Holder continuous, and only use the fact 
that Aper is Holder continuous (to obtain e.g. Lemmas [T51 [HI and ITTl) . The numerical examples that we have 
considered in Section [3] satisfy assumptions (|M | - ([55)) (remark that assumption (IM)) is also satisfied in the numerical 
examples considered in e.g. [50]). 

Note that we have assumed Ap^r and Bper to be symmetric only for the sake of simplicity. The arguments used 
below carry over to the non-symmetric case up to slight modifications. 
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4.1 Error estimation 



To bound the error between the reference solution and the MsFEM solution us, we use in many instances that 
we work in a weakly stochastic homogenization setting. We first recall in Section [4. 1.11 some results specific to weakly 
stochastic homogenization. This setting also allows to state rates of convergence for the two-scale expansion of u^, as 
we explain in Section 14.1.21 Our main result, Theorem 1 10[ is given in Section 14.1.31 



4.1.1 The homogenized equation 

Under the conditions recalled above, it is known (see e.g. [121 140]) that the solution u^(-,cj) to ((30)) a.s. converges 
weakly in Hq{'D) as £ to the deterministic solution u* of the homogenized equation 

- div [A*\7u*] = / in X>, u* = on dV. (36) 

The homogenized matrix is given by 



(37) 



where, for any p G S.'^, is the unique (up to the addition of a random constant) solution to the corrector problem 

-div [A^{-,uj){p + Vw;{-,uj))] =0 iuM^ 

VuiJJ is stationary in the sense of (j^ . , , 

(38) 



E 



^w;{y,-)dy^ =0. 



The variational problem associated with writes: find u* £ Hq{'D) such that 
where 

A*{u,v)^ I {\/v{x)f A*\/u{x)dx and b{v) = f f{x)v{x)dx. (39) 



•D 



As shown in [131 124] , in the weakly stochastic setting, the homogenized matrix A* can be expanded in terms of a 
series in powers of rj: 

A*^^A;^^+T^A\+7fAl{T^), (40) 
where ^2(77) is a deterministic matrix, that depends on 77 and is bounded as 77 — > 0, and where, for any 1 < «, j < d. 



= / (e,+VO^Ape.(e,+Vu;°^.), (41) 



'Q 

(Alh = /(e.+VO^E(Ai)(e,+VO, (42) 

where, for any p G R'', Wp is the unique (up to the addition of a constant) solution to the deterministic corrector 
problem associated to the periodic matrix Ap^r- 

-div [Aperip + Vw°)] = 0, 
Wp is Q-periodic. 

Under the assumption we have A* = E(Xo) B, with 

Vl<i,j<d, B,j^ [ {e, + Vwl fBper{ej+S/w°.). (44) 
JQ ' ' 

Remark 6. In general, when Ap^r is not symmetric, the expression of A* includes additional terms. Indeed, writing 

we in general nee(iE(V?«p) to compute Ai (see e.g. 124[\1^ ). In the symmetric case, these 
additional terms vanish, see e.g. JJ, Remark 4. 2 p. 117]. In the non-symmetric case, the expression (|42p of A\ needs to 
he slightly modified, but the expansion (j40p remains true. Our arguments hence carry over to the non- symmetric case. 



19 



Using the expansion (PU]) of A* with respect to 77, it is easy to see that the solution u* to ([5^ can also be expanded 
in a series in powers of 77. We have 



(45) 



(46) 



u* = u^ + T]E{Xo)ul + ifr^ with |!r^|!ffi(x)) < C, 
where C is a constant independent of 77, and where Uq solves 

- div [Al^^Vul] = / in 2?, = on 92?, 
and tl^ solves 

- div [Al^^Vu{] = div \BVul] in V, u\ = Q on c*!?. (47) 

The expansion (1451) will be useful in the sequel. We will also need a bound on u* and in the iJ^ norm. Recall that 
u* is the solution to (j36p . whereas is solution to 

- div [A*Vr^] = div [^^(77) (VwJ + 77E(Xo)V?Ii) + ¥.{Xq)A\Vu{] in = on dV. (48) 

In view of (1311) . we have, almost surely and almost everywhere, a_ Id < < a.)- Id in the sense of symmetric matrices. 
Recalling that homogenization preserves the order of symmetric matrices (see e.g. [48l page 12]), we deduce that 

In addition, the right-hand sides of (I36p and (1481) are bounded uniformly in 77 in the norm. Using |34l Theorems 9.15 
and 9.14], we obtain that there exists C such that 



Vr/, \\u*^\\m{p) < C and ||r^||i/2(x,) < C. 



(49) 



4.1.2 Two scale expansion of the reference solution n't 



As recalled above, the standard error analysis for the MsFEM in the deterministic setting is performed in the case 
when the matrix A"^ in ([Ij reads A^{x) = Apf,r{x/e) for a fixed periodic matrix Api^r- The problem ^ then admits 

a homogenized limit. To obtain bounds on the MsFEM error, one step of the proof is to approximate the oscillatory 

d 

solution by its two-scale expansion u* + diU* , where u* is the homogenized solution, is the periodic 

i—l 

corrector associated to p € M'^, and diU* — — — . In the deterministic case, it is known (see e.g. [T2j [23l |40] ) that, under 

OXi 

some regularity assumptions on Appr and u*. 



m{T>) 



for a constant C independent of e. 

In the stochastic case, it is known that 



a 

I* + Wei (~''^) 



(50) 



converges to as e — ?> 



(see [471 Theorem 3]), but no rate of convergence is known (except in some one-dimensional situations, see e.g. [101 1161 
143)). However, in the present article, and as announced above, we consider a weakly stochastic case. In this setting, we 
have derived in [32] a result similar to (ISUl) . We now state this result, which will be useful for our analysis. 



Theorem 7 (from [42] . Theorem 2). Assume d> \. Let he the solution to (1301) . and assume that A^, satisfies 
(ISS])-(I331)-(I3S]). Let A*^^, w°, Uo andll*^ be defined by (gt]), (|331), & and (gT]). The two-scale expansion of u"^ 
reads 

d 

v^^{-,uj) = u* + 7]E{Xo)ul + eJ2 [< (-) + vHXo)d,u'i) 

i=l 

+r,E{Xo)^e, (-) d,ul + 7jY, (^fe(w) - E(^o)) Xe, (- - d,ul , (51) 
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wher 



L = {fc e Z'^ such that e{Q + k) nV ^ ?)} , 



and where, for any p G M'', V'p ^■s the solution (unique up to the addition of a constant) to 



-divlAper'^ipp] = div [Bper {p + Vu;°)] 
tl^p is Q-periodic, 



(52) 



and Xp is the unique solution to 



-div[Aper'^Xp] = dw[lQBper(p + VwJJ)] 



Xp e Ll^{^% Vxp e {L^{W')Y' 

lim Xp{x) = 0. 
We assume that u* G W'^-°°{V) and u*^ G W'^^°°{V). Then 



(53) 



E 



(54) 



where C is a constant independent of e and rj. 



As pointed out above, and in [32], the assumptions p4p - ([55)) are standard assumptions when proving convergence 
rates of two-scale expansions (see e.g. [301 P- 28]). Likewise, the assumption Uq G W'^'°°{'D) (and subsequently G 
W'^'°°(V)) is a standard assumption (see e.g. Theorem 2.1] and |3D], p. 28]). In view of (j46p , this assumption implies 
that the right hand side / in (00]) belongs to L°°{V). 

In dimension d = 1, the boundary conditions of (j53p need to be modified for this problem to have a solution. We 
have derived in [32] the following result, which is the one-dimensional version of Theorem [7] (note that we need below 
weaker assumptions than in Theorem [7] as pointed out in [31]: we do not need to assume ((M]) -(|55 ]) . and the assumption 
/ G L'^iV) is enough): 

Theorem 8 (from |42j . Theorem 3). Assume that the dimension d is equal to one. Let uf^ be the solution to pop in the 
domain T> with f G L^(T>), and assume that A^j satisfies (j32p - (j33p . Let wf^ be defined by (I5ip . where the definition (j53l) 
is replaced by 

- [AperX']' - [l(0,l)Sper(l + (w°)')]' 

XGLLW, x'eL2(R), 



(55) 



where solves (j43p . T/ien 



E 



E 



(56) 



where C is a constant independent of e and rj. 



The following estimate, which is proved in [311 proof of Proposition 11] and useful to demonstrate ((M]) . will also be 
useful here: 

Lemma 9 (from W3 , proof of Proposition 11). We assume and d > 1. For any p G M'', any k £ , and any 
bounded domain T) C R'^, the function Xp defined by (|53l) satisfies 



Xp 



(--k) ' <Ce''Rd.e 



2 

for a constant C independent of k and s, where Rd,s = 1 if d> 2, and Rd,e = 1 + ln(l/e) if d = 2. 



(57) 
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4.1.3 Main result 

Before presenting our main result, we need some useful notation. Following the approach presented in Section [ 
recall that 

Wh ■■= span((?!)f , i = !,■■■ , L), 

where are the highly oscillatory MsFEM basis functions. By construction, the solution us G Wh of the weak 
stochastic MsFEM approach ((TC)) satisfies 

yvh€Wh, A'^^,^{us,Vh) = b{vh) a.s. (58) 

where, for any u and v in Wh, 

A^^{u,v)^ V / {\7v{x)f A.J-,uj)\/u{x)dx and b{v) = [ f{x)v{x)dx. (59) 

For future use, we also define, on the standard finite element space 

Vh :== span((?!)^, i = !,■■■ ,L), 

the forms 

Al^{u,v)^Y. I (V(7^k(«))(x))^A,(-,^) V(7^f,(^^))(a;)dx and Mv) = j f{x)n\^{v){x) dx, (60) 

where the local, linear operators TZ^ are defined on Vh by 

Vl<z<i, 7^k(0?|K)- '^flK- (61) 
These local operators give rise to the global operator TZ^ : Vh Wh defined by 

VK, yveVh, n%v)\^ = ni^iv\j^). (62) 

As pointed out above, the space Wh is not a subspace of Hl{T>), as the basis functions 0f may have jumps at the 
finite element boundaries (due to the use of the oversampling technique). We will hence work with the broken iJ^-norm 
introduced in (ITSl) . that reads, we recall, 

1/2 

'ivh^Wh, hhWrn " 



.KeTh 



We are now in position to present the main result of this article. We introduce the notation Qf = e{i + Q) for any 
i G lA, and denote by A'k the number of cells Ql in the element K: A'k = Card(i; Ql C K). We make in the theorem 
below a regularity hypothese on the macroscopic mesh, assuming that the volume of each element is bounded from 
below by ah'^, for some a > 0, and hence that TVk > a {h/e)'^. 

Theorem 10. Assume that Ajj satisfies (l32|) - (p3)) - ((34)) - ((35|l . We assume that Uq andu\ respectively defined by (j46l) 



and (|47p satisfy Uq G W^'°°{'D) and u* e W'^'°^{V). Let ufj be the solution to dSD]) and Us be the weakly stochastic 
MsFEM solution to (j58p . Suppose that d > I, e < h, and that there exists a > 0, independent o/K, h and e, such that 

Nyi >ai—\ . We then have 



^¥.\^\u-^~us\\l^]<c(^V~e + h+'-+ (1)'^' HN{h)) + 7/ + rfC{Ti)^ , (63) 

where C is a constant independent of e, h and rj, N{h) is the number of elements K in the domain T) (which is of 
order h^'^ in dimension d), and C is a bounded function as rj goes to 0. 

The restriction to d > 1 comes from the fact that the proof of this result uses the rate of convergence on the 
two-scale expansion of that we stated in Theorem [T] This rate of convergence is not optimal in dimension one, 
as can be seen from the comparison of (|5l)l and (j56p . The one-dimensional version of the above result is stated in 
Section [4.31 below (see Theorem I18L where we briefiy consider the one-dimensional situation. 

Remark 11. In the case rj = 0, our approach reduces to the standard deterministic MsFEM and we obtain the same 
estimate as in the deterministic case with oversampling (see e.g. \3'A Theorem 3.1]). 
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4.2 Proof of Theorem [JO] 

The proof of Theorem [TU] is the direct consequence of three lemmas. First we recall the second Strang's lemma 
Theorem 4.2.2, p. 210]. 

Lemma 12. Consider a family of Hilbert spaces Wh with the norm \\ ■ ||^i , a family of continuous bilinear forms ^ 
on W/i that are uniformly Wh-elliptic, and a continuous linear form b on Wh- For any h> 0, introduce us solution to 

yvh e Wh, A'i\^{us, vh) = b{vh) 

and G Hq{'D) solution to 

V«ei/i(2?), Al,{u-,^,v) = b{v). 

Then there exists a constant C independent of rj, h and e such that 

\\u^-us\\hI^<c[ \\uI^~vu\\hI+ sup ' ^'"^ ^' r ^ . (64) 



The first term in the right hand side of (|64p is the so-called best approximation error. The main part (step 2) of 
the proof of Theorem [TU] is devoted to its estimation, following up on the estimate ((M)) provided by Theorem [T) 

The second term in the right hand side of (1641) is the so-called nonconforming error, which vanishes in the case 
Wh C Hq{'D) (the method is then conforming, and we are left with the standard Cea lemma). In our case, we use the 
oversampling technique, hence our approximation is not conforming, and this second term does not vanish. It will be 
estimated in the step 3 of the proof of Theorem 1101 using the following two results, which are proved in Appendix \K\ 

Lemma 13. Consider the two bilinear forms A* and A^ respectively defined in i39\) and I160\). Under assumption p4p . 
there exists a deterministic constant C, independent of rj, e and h, such that, for any Vh & Vh, 



sup 



A^^jj{vh,Wh) - A*{vh,Wh) 



< C 



ill ^ ^'^^ ^ 'n^^i'n)) hh\\m{v) a.s., 



where C is a deterministic function independent of e and h and bounded when rj 



h, e) 



max max 

K l<p,m<d 



1 



:, + V<(i)f(.4,(i.)-E(A,( 



where is the largest domain composed of cells of size e included in K; 



0, and A is defined by 

t (?) 



dx 



(65) 



(66) 



rCK 



Lemma 14. Consider the two linear forms b and bh respectively defined in fi39[) and Ii6&[). Under assumption p4p . 
there exists a deterministic constant C independent of rj, e and h such that 



sup 



bh[wh) - K'^h) 



\\wh\\m{v) 



<Ce\\f\\L^vy 



(67) 



Before turning to the proof of Theorem [TUl we first give some properties of the random variable X{uj,h,e) that 
appears in the right hand side of (|55|) , and we next detail a two scale expansion of the highly oscillatory basis functions 
which will be useful in the sequel. 

Remark 15. We will show in Lemma \TE\ below that X defined in (j66p is uniformly bounded with respect to h, e and 
LO . Since A* is coercive, we deduce from (|65p that A^ ^ is also coercive, in the S( 
constant a > 0, independent of h, e and rj, such that 



that there exists a deterministic 



Vu/i e Vh, a\\vh\\H^xi) < A'^^^{vh,Vh). 
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4.2.1 Properties of A(w, /i, e) 

We state here some useful properties of the random variable X{uj,h,s) that appears in (|65|) . They will be proved in 
Appendix [21 As mentioned above, we recall that the assumption Nk > a. (h/e)'^ that we make below is a regularity 
assumption on the macroscopic mesh (the volume of each element K is bounded from below by ah'^). 

Lemma 16. Let \{u)^h,e) he defined by (j66p . Then, there exists a deterministic constant C such that, for any h and 
£, we have < A(w, h,e) < C almost surely. 

Assume now that the random matrix Ai satisfies (j33p . where the law of Xk{uj) is absolutely continuous with respect 
to the Lebesgue measure. Assume furthermore that the number = Card(i; Ql C K) of cells in K satisfies -/Vk > 



a y~ j f^''' some a > independent of the element K, h and e. Then 

E{X{;h,ef)<C^[HN{h))]\ (68) 

where, we recall, N{h) is the number of elements K in the domain T> (which is of order h^'^ in dimension d) and C is 
a deterministic constant independent of h and e. 

Because of the specific form of A\ , we will see in the proof of that result (see Appendix [2] below) that 

A(a;, /i, e) = max max , (69) 

K l<rn,p<d 

where each random variable S^'^ is a normalized sum of (h/e)'^ i.i.d. variables. Applying the Central Limit Theorem, 
we hence know that 5*^'^ converges, when e — > 0, to a Gaussian random variable (up to an appropriate renormalization). 
Likewise, computing the expectation of (5'^'^)^ is not difficult. However, in the above lemma, the difficulty stems from 
the fact that A(cij,/i,e) is the maximum of many such random variables S*^'^ (in (j69p . the number of elements K is 
indeed of the order of h^"^). Our main task is hence to control how E(A(-, h, e)^) depends on h. See also Remark \W\ 

4.2.2 Two-scale expansion of the highly oscillatory basis functions 



Following [32], we recall here an expansion of 
we have, for any 1 < i < d + 1, 



where an is such that 



that will be useful in the sequel. By definition (see PT|) and (IT^ ). 

d+l 

6,S 

K 



(70) 



iO,K 



d+l 

i=i 



o,s 

j 



K 



We thus first turn to Xi'^ ^ which, by definition (see (1101) ). is the unique solution to 



div 



Aper (J) Vxr^(a;)] = in S, y^^^^{x) - xl'^ix) on 9S. 



We introduce the function 



(71) 



(72) 



(73) 



where w". is solution to the periodic corrector problem (|43|) . By construction, using (l73t . (|43)) . (j72p and the fact that 



Vx^'^ is constant on S, we have 



—div 



A 



per 



(?) 



^ef{x) 



= —div 



-div 



£ 

0, 



V l^xf (x)-x°''(x)-s5:< Qa,x° 

Ar>er (f ) Vxr^(x)] - i^a.x^'div [a,,. (^) (e, + V< ( 



i=i 
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a 

while, from ([75]) . 6l'^{x) = djXi'^ {x)'w^. (^—^ on dS. So, by linearity, we obtain 

d 

^■'^(^)=E^^-^^'^'(^)' (74) 

where g ^^^(S) is the unique solution to 

- div [Aper (j) V<e^(x)] =0 in S, = (J) on OS. (75) 

Using ([751) . now obtain a useful relation between (j)^'^ and Indeed, collecting (ffn)) . ([7T|) . ((75| and ((71)) . we 

obtain the exact expression 

<^r^(-) = 0^ (-) + (<. ) - ^i{x)Q 9,0^- (76) 

Recall now that 0^'^(x) = '^k('^i'^)j by definition of the local operator TZj^ (see ((6T|) ). Correspondingly, the global 
operator TZ"^ , defined on Vh by (p^ , equivalently writes 

d 

W € V,, 7e^(u) = + (< (-) - ee) d,u, (77) 

where is locally defined on each element K by = ^eIk' construction, for each K, £^1 e i/^(K), but it a 

priori does not belong to H^{'D). The relation ((77)) allows to extend the operator TZ'^ on H^(T>). 

We now recall the following bound on the function , that appears in (j76p . In [32], this lemma is stated in dimension 
d = 2, but its proof, which essentially makes use of [9j Lemma 16], carries over to any dimension. 

Lemma 17 (see 05], Lemma 2.1). Let he the solution to ([75]) . with Aper satisfying ([M)) . Consider K C S, with 
diam(K) — h and dist(K,dS) > h. Then there exists a constant C independent of h and e such that 

l|VeilU~(K) < J. (78) 

4.2.3 Proof of Theorem [lO] 

The proof is based on the bound (|64p in Lemma [T^ where the bilinear form ^ and the linear form h are defined 
by (|59]) . In Step 1, we show that the bilinear form is coercive for the norm \\ ■ defined by (fT8)) . Step 2 is 
devoted to appropriately selecting an element Vh £ such that — I'/iHif i can be analytically estimated. This will 
provide a bound on the first term in the right hand side of (|64p . In Step 3, we bound from above the second term in 
the right hand side of ([M)) . using Lemmas [T51 and [TH Step 4 collects our estimates and concludes. 

Step 1: We first show that the bilinear form A^ ^-^ defined by (15^ is coercive for the norm || • ||^i defined by ([T5)) . 

Consider the bilinear form A'^ ^ defined by (|60p. We pointed out above (see Remark [T5| that it is coercive on Vh- 
Hence, there exists a > such that, for all Vh G Wh, 

a\\vh\\m(v) < ^lr,{^h,vh) = A!l^{vh,vh), (79) 

where Vh € Vh is such that Vh — TZ^(vh)- Since, in the bilinear form A'^ ^j, the matrix is bounded, we deduce the 
estimate 

WvhWl^v) < C\\vh\\l., (80) 

that we will use in the sequel. The sequel of this step is devoted to proving that there exists C independent of h and 
e such that, for all Vh E Vh, 

WvhWni < C\\vh\\Hi{v) ^^thuh^Wivh). (81) 



25 



Combined with ([7^ . this shows that A'^ is coercive for the norm || • 

To prove ([ST|) . we first write that, since Vh — TZ^(vh) and Vh € Vh, there exist some coefficients {Pi}^^^ such that, 

for any x ^V, Vh = X^i^i A'/'i ^'^d ''^h = Ti-^{vh) — Si^ift'/'f- Consider now an element K, and its corresponding 
oversamphng domain S. We know from ((TT]) and (|13p that 

L d+l L d+l 

i=l j = l 1=1 j=l 

Consider now the functions 

L d+l L d+l 

i=l j = l 2=1 J = l 

defined on S, and that satisfy, by construction, 

VxeK, Vh{x)^w^{x), Vh{x)^wl{x). (82) 

In view of ([TU)) . we have 

-div [A^(a;)Vu;f (a;)] = in S, wf=wf on dS, 

which implies that 

\\wf^\HHS)<C\\wt\\HHS)- 

We deduce from (|82p and the above bound that 

ll«ft||Hi(K) = llw'fllHMK) < I|w^||h1(S) < C\\wf\\H^S)- (83) 

We next see that there exists C independent of h such that, for any piecewise-affine function r on S, we have ||T||/fi(s) 5: 
C'|lT||/fi(K)i provided there exists < c_ < c+ independent of the element such that c_ < < c+. Using this bound 
for T = wf, we infer from ([83]) and ([82]) that 

||wh||Hi(K) < C||wf ||h1(s) < Cllwf ||_ffi(K) = C'|lw/i||_ffi(K)- 

Summing over all elements K, we obtain (j8ip . and this concludes this first step. 

Step 2: Let U'^u* be the projection of u*, solution to (1551) . on the standard FEM space V/j. We have TZ"^ (ll'^u*) S 
W?i (recall TZ'^ is defined by (|62|). and equivalently writes as in ((77)) ). Our argument is based on the following triangle 
inequality: 

^{^%H~-'^\\k) ^ ^^{H--rrHHv,)+2^[^%H-v&^ (84) 

(iK, - <Xhhv)) + 2IE (ll^^ - 7^^(^^i;)||2^,) (85) 



< 2 

< 2: 



■E (h^, - v'r,\\m(v)) + 4E (^^ - 7e^«)||?^i) + 4||7^^«) - 7^^(^''^.;)||^, , (86) 

where v^{-,Ld) G H^{'D) is defined by (|5ip . The estimate ([5^ in Theorem [7] bounds the first term from above. In the 
following two sub-steps, we bound the other two terms of ((55)) . 

Step 2a: bound on E - 7^^(^i*)||^J 

Using the expansion (P5|) of u* in a series in powers of 77, and ([77]) . we write 

d 

n^{u*) = u* + riE{X^)u\ + e ^ (-) ~ If) (a^uj + 7?E(Xo)aX) + ^'5,, 

p=l 
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where 



Using (l5lT) . we thus have 



p=i 



(87) 



z;^(-,c.) - 7^^«,) = 77£^ ( E(Xo)^.e, (-) Sp^/S + ^ (Xfc(c.) - E(Xo)) 



Xep ( - - ^ ) dpul 



eJ2&idpu;~v'dpr,)^Tj'g,. (88) 
p=i 



To bound E 



K 



, we first establish a few simple results. First, there exists 5 > such that, for 



any 1 < p < d, we have 



(89) 



where Wp is the periodic corrector, solution to (|43|) . This is a consequence of the fact that Ap^r is Holder-continuous 
(see assumption ([M])). in view of [34j Theorem 8.22 and Corollary 8.36]. We infer from (|89p and the periodicity of w'^ 
that, for any 1 < p < d, we have 

^0 g ^l.oo(]R'i). (90) 



Second, for any 1 < p < d, we have 



Xe 



, (- - fc) ' < Ce''Rd,e and Vxe, e (i'(K'^))'' , 



(91) 



where C is independent from e and fc, i^^.j = 1 if d > 2 and Rd,e = 1 + ln(l/e) if d = 2 (see ([55)1 and ([ST]) ). Third, we 
see that, for any I < p < d, 

C 



<C and llVeiLoo(K) < ^, 



(92) 



where C is independent from e and h. The second assertion is given by Lemma 1 171 above, whereas the first assertion 
comes (ESI): using again [Ml Theorem 8.22 and Corohary 8.36] and we first see that, for any S, € C^'^{S) 



for some 6 > 0. Using next the maximum principle on (j75p . we have ||C?IIl°°(s) ^ f^cp ^-^ 
using dH?]), dUl), dMl) and ([Ml), we obtain that, for any element K, 



< C. Lastly, 



I|5i)IIhi(k) < C'lk»?llH2(K) (1 + ^ + /jj - c'II^»jIIh2(k), 



hence 



E 



hvWh^-o) < XI II5'?IIhi(k) ^ C!\K\\h2(^v} < C*. (93) 

K 

We are now in position to estimate (gll). Using that G W^'°°{'D), we deduce from and dM]) that 

H-nHu;)\\l.^J < Cr^^e'Y.\\dpU*o\\U (E{Xor ^eJ-) ' ^ + Var(Xo)5]||xep(--fc) ]) 

p=i \ ^ fce/, ^ ^ ^^V 



< C7]^e^ (1 + (Card 4) £'*i?d,e) + Ce^ + Ct;^ 



(94) 



for some constant C independent of e, t] and /i, and where, we recall, i?^,^ = 1 + ln(l/e) if d = 2 and Rd,e = 1 if d > 2. 
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We thus have a bound on u^(-,a;) — Tif(u^) in the L"^ norm. To prove a bound m the broken norm, we consider 
V?;^(-,cj) — WTZ'^{u*): we see from that, in each element K, 



(95) 



where 



Dn = 



D2 = e 



HWe, (-) ^dpUl + ^ (Xfe(c^) - E{Xo))Xe 

mxoW^e, (-) + J2 (^fe(^) - E(Xo))vxe, (- - fc) apw5 



Note that Z?o and Di are globally defined on V, but _D2 is not (as ^1" may have jumps from one element K to the 
other). We now bound these three quantities. Using and the fact that Uq G W'^'°°{'D), we have 



I 2 

E(|Polli.(p)) < CYW^dpKWl- (nXo)' ^ej-) +Var(Xo)^ 



(96) 



where C is a constant independent of e and h. We now turn to Di: using again (PTt and that G W^^'°°(2?), we 
obtain 



l^il 



p=l \ 
< (l+Var(Xo)^e'^||Vx, 



2 



-Var(Xo)^ ||vxe, (--A:) 



5=1 \ fce/s: 

where C is a constant independent of e and h. Turning to D2, using (j92p . we have in each element K that 



(97) 



hence, using P^ . 



||i?2||i2(K) < Ce'K - ?]V,||lf2(K) (1 + /J2) ^ - 

2 2 



2„ l|2 

//2(K)J 



(98) 



Collecting ([M]), (IMl), (ISZ]), (IMl) and ([Ml), we obtain that 



E 



El|v<,-v7^^(u;)|li.(K) 



K 



CoUecting this bound with (IM|) . and assuming that I?]] < 1 and e^Rd.e < 1, we deduce that 



E 



|^;^~7^-(^;)||^. 



= E 



K 



(99) 



where C is independent from e, h and 77. 
Step 2b: bound on \\n%u*) - W (Il''u*)\\jj^ 
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Recall that TL^u* is the projection of u* (on the standard Pi FEM space Vh), hence |jn''u*|| jyicp) < \\u*\\h^{v)- 
In addition, using l^^ . we have, using a standard result from the theory of Pi finite elements (see [2H Theorem 3.1.6 
p. 124]) 

IK, - n"u;iU2(P) + h\\u; - n'^u^Wmiv) < ch^s/^u;\\mv} < ch\ (loo) 

where C is a constant independent of h and 77. In view of (|77p . we have 

d 

n'{u^^) - 7e^(n''7.;) - < - u^u*^ + (< (-) - ~ n"u;). 

p=i 

We deduce from dHO]), (El and (ITU(I|) that 

||7^^(w;) - 7e^(n''7.;)|U2(p) < c{h^ + e/i). (loi) 



We now turn to bounding the gradients. Recall that V(n''M*) is constant in each element K. We thus have, using ([90 
and JHl), that 



|V7^^(«;) - v7^^(^''^i;)|U.(K) < IK - n''<,lkMK) i + E ( ||^< (-) 

\ p=i 



L=(K) 



-£||veilL-(K; 



t/i^,+ I 



L~(K) 



,,ll-ff2(K) 



< c|K, - n'^u;i|Hi(K) (1 + + f iisii^^^(K)- 

We then deduce, using (ITUg|) and dH]), that 

^ ||V7^^(M;) - v7^^(^''u;)||i.(K) < cn^; - n'^u^^^.^T,) + + ^'IK'l^^c^') 

K 

< C/l^ (l + ^)'+C£2, 

where C is a constant independent of e, 77 and /i. Collecting this bound and (jlOip . we obtain 



(102) 



K 



where C is a constant independent of e, rj and h. 

Step 2c: We are now in position to bound the first term in ([M]). We infer from ([HS]), dMl), (|M1) and (|102p that 



E inf ||ii^-w/.||^i < C + ?7v/e ln(l/e) + + - + /i 



(103) 



where we have assumed that eln(l/e) < 1. 



Step 3: We next turn to estimating the non-conforming error, namely the second term of the right-hand side of (|64p . 
For any Wh G W/i, introduce Wh & Vt such that TZ^{wh) = wt (recaU that TZ"^ is defined by ((62)) ). We note that b{'Wh) = 
bhiwh), where the linear forms b and bh are defined by (j59p and (1601) . Using the weak form of the homogenized equation 
(see ([M]) ). we see that b{wh) = A*{u*,Wh)- In addition, by definition of A^^^ (see (pO])). we have A^.^(nhU*,Wh) = 
A^,^(n''{UhU*),Wh). For any wh G W/i, we have 

\A^^^{u'^,Wh) - biwh)\ < \A'^^^iu'^,Wh) - A'l^in'{llhu*),wh)\ + \A'^^^in'{Ilhu;),Wh) ~ b{wh)\ + \b{wh) ~ b{wh)\ 

< ||A^||loo II - 7^''(^,i^t*) 11^:^1 \\wh\\Hl + ■^'lr,(^hU*,Wh) - A*{u*,Wh) + b{wh) ~ bh(wh) 

< \\AjL^\\ut^~n'{iihu;)\\Hi\\wh\\H}^ + Ai\{ii''u;,wh) - A:^{ii''u;,wh) 

+ \\u*-U''u*\\Hi(v)\\wh\\HHv)+ b{wh)-bh{wh) 
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where we have successively used the continuity of the bihnear forms ^ and . Using Lemmas [T3] and [U for the 
second and the fourth terms respectively, we deduce that 

\Ai^{u'^, Wh) - b{wh)\ < C\\u'^ - n'{IlhU*)\\Hi\\wh\\Hi + ^ (J + VM.^^ h,e) + TfC{ri)^ \\n''u*^\\mip)\\wh\\m{T>) 

+C|K, - ^''K,\\HHV)\\wh\\H^(v) + Ce\\wh\\H^[v), 

hence, using ([501) and (ITUg| . 

\Al^{u-^,WH)-h{wn)\ ^^„^^, . , , = , , . .^ , . , .2^..^^ (104) 



<C|l<,-7^^(^,^i;)|l^l+C /i+-+7?A(c.,/i,e)+e + 7?^C(7y) 
The first term is bounded as in Step 2; 

IK, - 7^^(^,<,)||^l < IK, - v^^Wh^v) + IK, - t^'{K)\\hi + WiK) - 7^^(^,<,)||^l 

hence, using (|54l) . (IMl) and (I102p . and assuming that eln(l/e) < 1, we have 



E 



Collecting (|104p and (jlOSp . we thus obtain 



E 



sup 

.wheWh 



<C{h+-+ 77\/E [A2(-, h, e)] + + r; + r?2c(,y) 



(105) 



(106) 



Step 4: Collecting ([Ml), PII5|) and PIM)) . we get 



E [||Mf, - usll^i] <C7(Vi+/i+^ + r/VE]A2(TM] + ^ + '7'C(77)) 

where C is a constant independent of e, /i and ?/, and C is a bounded function as 77 goes to 0. Using (|68|) . we deduce 
that 

^eJiIu^- 7.511^1] < C (^V^ + /I + ^ + (i) ln(7V(/.)) + + 7?2c(r;)^ 

where A''(ft.) is the number of elements K in the domain (which is of order h^'^ in dimension d). This concludes the 
proof of Theorem [TOl 

4.3 The one dimensional case 

In this section, we briefly consider the one dimensional situation. As in the multi-dimensional case, we assume here 
that afj(a;, w) — a,, ^— , , where a,, is a stationary random function satisfying, for any I77I < 1, the condition < a_ < 
0^(2:, w) < almost everywhere in R, almost surely. In line with p2p. we assume that 



a,,(x,a;) = aper{x) + 77 ai{x, cj). 



(107) 



where 77 is a small parameter (I77I < 1), ap^r is a 1-periodic function satisfying the condition < a_ < aper{x) < 
almost everywhere on M, and ai is a bounded stationary random function: |ai(a:;,w)| < C almost everywhere in M 
almost surely. In the vein of we suppose that 



ai 



{x,u!) — ''^l(^k^k+i]{x)Xk{oj)bper{x) such that 3C, Vfc G Z, \Xk{oj)\ < C almost surely, (108) 



where {Xk{uj))i^^^ is a sequence of i.i.d. scalar random variables and bper G L°°(K) is a 1-periodic function. 
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Note that, in this one-dimensional setting, we do not make any regularity assumption on Uper (in the vein of ([M])). 
In the multi-dimensional case, this assumption is useful to e.g. state that the periodic corrector satisfies Wp G W^'°°{M.'^) 
for any p S R''. In the one-dimensional case, the corrector problem can be solved analytically, and one can see that the 
above assumption aperix) > a_ > almost everywhere on M is sufficient to obtain such regularity on the corrector. 
Similarly, we do not need to assume here, in contrast to Theorem [TUl that Uq and u* defined by (pS)) and (|47l) both 
belong to W'^'°°{'D) (an assumption equivalent to / G L°°{'D), in the present one-dimensional setting). The assumption 
/ G L'^{V) is sufficient. 

The problem (l30l) now reads 

-^(^a^[^,u:)^ut,{x,uj)^^f{x) in (0,1), ^(0, ^) = w^(l, t^) = 0. (109) 

We consider a uniform discretization of the interval (0, 1) in the elements = {xi,Xi^i), with x^+i — Xi ~ h = 1/ L 
for some L G N*. 

The one-dimensional version of Theorem 1101 reads as follows: 

Theorem 18. In the one- dimensional setting, assume that af^ satisfies (jl07p - (|108p . Let u,^ be the solution to (I109P 
with f G 2.2(0, 1), and Us be the weakly stochastic MsFEM solution to (j58p . Suppose that h/e G N*. We then have 



yE[hf,-us||^,] <c(^e + h + rj(^^) ^'^ \^{\lh) + + ifC{ri)\^ , (110) 
where C is a constant independent of e, h and i] and C is a bounded function as rj goes to 0. 

Proof. The proof of this result follows the same lines as that for the multi-dimensional case. It is based upon the 

homogenization result contained in Theorem [5] above. □ 

As pointed out above, the rate of convergence stated in Theorem[7] (and hence the estimate provided by Theorem [TU|) 
is not optimal in dimension one. This hence motivates Theorems [5] and [T51 which are their respective one-dimensional 
variants. On another note, the assumption h/e G N* implies that some terms in the error bound vanish. A result 
similar to (IllOp holds in the absence of such assumption, with the additional term e/h in the right-hand side. 



A Proofs of Lemmas I13L [14] and [16 



Proof of Lemma \13[ This result relies on the expansion 

d 



^(.) = 0°'^(x) + e y: (f ) - er(-)iK) 



m— 1 



from (|76p and the fact that V^^' is constant on K. 
For any Vh and Wh in Vh, we write 



Y if iVwHix)fA;VvHix)dx~Y«f (v#^)^A,(^,a.)v</.^^-^da; 
Ken hj=i •'^ J 



where Vh — ^ ^h't''j likewise for Wh- Using the above expansion of cjf, ^ and the fact that Vcj)"''^ is constant on 



iO,K ■ 
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K, we have 



d 



E 

m,p—l 



1 

|K| 



'K 



K 



= E 

m,p— 1 

We thus obtain 



1 

|K| 



K 



da; 



- £N 



K 



E E / ^rnVhdpWh 



< EI1^''II^MK)I|w/i|Ih1(K) E i^™pi 



(111) 



K 



where 



mp |K| 



K 



dx. 



which we write 
with 



= Do + Di- D2, 



D2 = 



IKI 



IKI 



K 



(112) 

(113) 
da; I , 



K 



(ve(x))^A,(pw) ve(a;)dx. 



We are thus left with bounding |A^p| from above. We first bound Di and Z?2- Using Lemma [T71 (recah that Ap^r 
satisfies ([M)) . i.e. is Holder continuous) and the fact that G H^{Q) and is Q-periodic, we obtain 



and 



hence 



e 



K 



ep + V< Ver(a;)dx 



< 



< C 



< c£ 



\K\h 

e 
h 



(^|K|+e'^j|^^jV<(y)|dy 



(114) 



(115) 



where C is a deterministic constant independent of h, e and t]. We next turn to Dq. We introduce the cells Qf 
e{Q + ie U^, let 11^^ [j Qf, and recast pi^ as 

QfCK 



(116) 
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with 



bulk 



^boundary ^ 



+ V<(f)]"A,(f,.) 



da;. 



We denote by the set of cells Qf that intersect the element K, i.e. 

Jk^ U Q'- 

By construction, '^k- Using that G H^{Q) and is Q-periodic, we write 

1 



boundary 




IKI 



da; 



< C 



e'* |aK| 



K| 



ep + VwO^ (y) dy W / [e,„ + Vi(;0„ (y)] ' dy 



< C- 



We next consider (jll7p : 



with 



bulk 


|K\/&| 





|K| 









mp IJ^ 



'-D 



' 



Using the expansion (j40p of A* , we write 

r 1 ^ 

ep + (y) ^per (y) [em + Vw°^ (y)] dy 



—bulk 

= 









\Ik\ 


L 


Gp 






ep 








\ik\ 



A 



per 



t (?) 



da; 



(117) 



(118) 
(119) 



=P + V< [e™ + V<„ (^)] dxj +v^C{rj). (120) 



The leading order term in (|120l) vanishes. We are hence left with 



—bulk _ ri 



ep + Vu;° 



E Ml 



Vw 



Collecting ([TT^ . (fTTi)) . (ITT5|) . (fTTC)) . (fTTSl) . (ITTO| and (|nT|) . together with the fact that 



i^{^)]dx + rj'C{^). (121) 



K 



obtain 



lA^pl < C 



=, + V<Q]^(E(..(i.))-..(i.)) 



[An] 



dx 



mp 



< C-, we 
h 



We set 



Xiuj.h.e) = max max 

KeTfc l<m,p<d 



Kl 



ep + Vu;:?^ (^) 



E(A,(f,.))-A,(f,.)) 



da; 



33 



a 

We thus have, for any K, ^ | A^^l < C + rj\(uj, h, e) + rj^Ciif)^ . Using (jllip . we thus obtain that 

m,p— 1 

< +?7A(a;,/i,e) +77^C(7?)) Ilwftlkvx') 

This concludes the proof of Lemma [TS] □ 

Proof of Lemma \lJi\ Again, as for Lemma I13[ this resuh rehes on the expansion (|7S)) of 0- ' and the fact that V4>1' 
is constant on K. 

Setting Wh{x) = ^w\(f)^{x), we observe that 



hh{wh) - b{wh) 



K 1=1 



Using (1751) and the fact that V(/>j' is constant on K, we obtain 
L „ d 



We have 



and 













/ dpWh 


^=l 







< V|K| ||w/i||ifi(K) 



(123) 



Jji^) (< (f) dx\ < ||/|U.(K) (-) 



L2(K) 



ieiU^(K; 



< 



||/||i.(K)/^(||<L=o(R.) + ||eiU=o(K)) . (124) 



Recall now that, since Aper satisfies ([M)) (i.e. is Holder continuous), we know that and are both continuous, 
and that G L°°(M'^). Using the maximum principle on ([75]) . we write 

lieili-(K) < ll<llL-(as) < II<IIl~(r^), 

and we thus deduce from (|124p that 

./(^) (<(f)-eK 

for a constant C independent of h and e. Collecting (|122p . (|123p and (|125p . we obtain 

6/i(w/i) - 6(w/i) < C£||/||L2(x,)||u;/i||_ffi(x,). 



K 



<q|/IU.(K)v^ 



(125) 



This concludes the proof of Lemma [Ml 



□ 



Proof of Lemma \16[ We first prove the uniform bound on A. Recall that the field Ai is bounded almost surely and 
almost everywhere. This implies that 



1 



'Kl Ji\ 



+ V< )]" {a, (f - E {a, .))) [e„ + (f ) 



<2Mi|Uo 



'Kl 



ep + V<(-) 



e™ + Vw°„, (- 
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Then, using the Q-periodicity of , we obtain 



e^ + Vwl (-) - "E / k + Vii;° (-) 



dx = \n 



Kl 



We thus have 



X{uj,h,e) < 2\\Ai\\l°-=- max 

l<p,m<d 



hence X(oj,h,e) is bounded almost surely by a deterministic constant independent of h and e. 
We next turn to (l68l). Rewrite (l66l) as 



A(a;, /i, e) = max max |S'^'''| , 

K l<m,p<d 



with 



Kl J I 



+ V<(f)f(^.(i,.)-E(A(i,.))[e„. + V<,(i) 



Using the periodicity of the correctors and the specific form ([55)) of Ai , we have 



- E(Xo) 
K.;Q^/^ VVar(Xo) 



(126) 



with 



T™^P = VVar(Xo) /_ ep + Vu;°^ (y) Bp.riv) [e™ + (y)] dy and TVk = Card{i; C /^}. 

Thus, A(cj,/i,£) reads 



A(a;, /i, e) = 7max |S'k(<^)| , 



K 



where 7 = max r™'^ and 

l<m,p<(i 



Introduce the probability density function ip n-^ of the random variable \^N-k\ Sf^ {<^)\, and -Fjvk (2;) = P (V-^kI -^k I ^ 2;) . 

Using the assumption that each element K contains a number A'k of cells of size e that satisfies A^k > a f — j for 
some a > 0, independent of K, /i and e, we write 



E a 



< 



E(7VKmax|S'^p) = /" x^-^P (v^max |S'|c| < dx. 



Since 



we deduce that 



(^y]VKmax|S'K| <x] = 



N(h) 



N{h) 



E a 



< / X^N{h)Fl:l'-''^'\x) ^N^{x)dx = 61+62 



where 



61= I x^N{h) F^l'^^~^{x) LpM-i^{x)dx and 62= / x^ N{h) F^^^^^^ {x) lpnk{x) dx, 



(127) 



with Ch = 2\n{N(h)). We now successively bound from above ei and 62- First, integrating by part, and using that 
< Fn^ < 1, we obtain 



< ei = 



X = Ch, 
£C = 



2xFZ\x)dx<cl 



(128) 
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Second, again using < Fn-^ < 1, we get 

< 62 < r x'N{h)^N^{x)dx = N{h)E (lr^^^^^.NK\Sk\' 
Using the Cauchy- Schwartz inequahty, we obtain 

el < N{hf E [Nl\SI^\^] P (0Vk|^kI > cn) ■ 

Introduce Yi = — ^^====r, so that S^iuj) = Yi{uj). Recall now that {Yi)^^^d is a sequence of independent 

identically distributed variables, with mean zero. Any such variables satisfy the bounds 



Vp e N*, 3Cp > 0, VA^ e N*, 



E 



2p 



< 



Cp 



for a constant Cp that depends on p and the moments of Yi, up to order 2p. Recall that all moments of Yi are well 
defined, as Yi is bounded almost surely. Thus 



ei < a 



(129) 



We now recall the Markov inequality: for any positive non-decreasing function ip on R, and any real-valued random 
variable Z , we have 

miz)) 



V6 e M, F{Z >b)< 



We apply this inequality to the random variable Z{uj) = VAk'S'k(w), with ip — exp(<:-) for some t > 0, and b = Ch.. 
This yields 



exp 



K 



(130) 



where we have used the fact that is a sum of i.i.d. variables. Using a Taylor expansion with respect to t, we see 
that 



E 



exp 



t 



--Yo 



Thus 



K 



exp 



1 



K 



K 



-mo' 



1 



6N, 



3/2 
K 



E 



r^* exp(^yo/ VAk) for some ^ G (0, t) . 



< exp 



|E(yo') + — ^E (yo3exp(eyo/v^) 

2 ovAk ^ 



Using (|130p . taking t = 1, and using that e ~ , we obtain 



< 



N{hy 



N{hy 



■ exp 



■ exp 



N{h)^' 

2 ^ °^ 6VAk 
^EiY,^) + h{\Yo\'eM\Yo\)) 



EiYo'eM^Yo/^NK] 



for some ^ G (0, 1), 



(131) 



Likewise, considering Z{uj) = — \/AK5'|j-(a;), we obtain a similar bound. Collecting (|129p . (|13ip and the fact that Iq is 
bounded almost surely, we have 

el < C, (132) 
with C independent of h and e. Collecting (jl27p . (jl28p and (jl32p . we get, for a constant C independent of h and e, 

Ei\i;h,er)<C^[HN{h))f. 
This concludes the proof of Lemma [121 D 
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Remark 19. The above proof shows that, when e — > 0, the random variable y~ j X{i^,h,e) converges in law to 
= max IGkI^^^)!, where G-k{uj) are i.i.d. Gaussian random variables. Precise results on the behavior of Qhii-^) 
when /i — > (i.e., when the number of Gaussian random variables involved diverges) are given in e.g. 14 6} - 
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